Intrepid
test_02.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
56#include "Teuchos_oblackholestream.hpp"
57#include "Teuchos_RCP.hpp"
58#include "Teuchos_GlobalMPISession.hpp"
59#include "Teuchos_SerialDenseMatrix.hpp"
60#include "Teuchos_SerialDenseVector.hpp"
61#include "Teuchos_LAPACK.hpp"
62
63using namespace std;
64using namespace Intrepid;
65
69
71void rhsFunc(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
72 if (degree == 0) {
73 for (int cell=0; cell<result.dimension(0); cell++) {
74 for (int pt=0; pt<result.dimension(1); pt++) {
75 result(cell,pt) = 1.0;
76 }
77 }
78 }
79 else if (degree == 1) {
80 for (int cell=0; cell<result.dimension(0); cell++) {
81 for (int pt=0; pt<result.dimension(1); pt++) {
82 result(cell,pt) = points(cell,pt,0);
83 }
84 }
85 }
86 else {
87 for (int cell=0; cell<result.dimension(0); cell++) {
88 for (int pt=0; pt<result.dimension(1); pt++) {
89 result(cell,pt) = pow(points(cell,pt,0), degree) - degree*(degree-1)*pow(points(cell,pt,0), degree-2);
90 }
91 }
92 }
93}
94
96void neumann(FieldContainer<double> & g_phi, const FieldContainer<double> & phi1, const FieldContainer<double> & phi2, int degree) {
97 double g_at_one, g_at_minus_one;
98 int num_fields;
99
100 if (degree == 0) {
101 g_at_one = 0.0;
102 g_at_minus_one = 0.0;
103 }
104 else {
105 g_at_one = degree*pow(1.0, degree-1);
106 g_at_minus_one = degree*pow(-1.0, degree-1);
107 }
108
109 num_fields = phi1.dimension(0);
110
111 for (int i=0; i<num_fields; i++) {
112 g_phi(0,i) = g_at_minus_one*phi1(i,0);
113 g_phi(1,i) = g_at_one*phi2(i,0);
114 }
115}
116
118void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
119 for (int cell=0; cell<result.dimension(0); cell++) {
120 for (int pt=0; pt<result.dimension(1); pt++) {
121 result(cell,pt) = pow(points(pt,0), degree);
122 }
123 }
124}
125
126
127
128
129int main(int argc, char *argv[]) {
130
131 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
132
133 // This little trick lets us print to std::cout only if
134 // a (dummy) command-line argument is provided.
135 int iprint = argc - 1;
136 Teuchos::RCP<std::ostream> outStream;
137 Teuchos::oblackholestream bhs; // outputs nothing
138 if (iprint > 0)
139 outStream = Teuchos::rcp(&std::cout, false);
140 else
141 outStream = Teuchos::rcp(&bhs, false);
142
143 // Save the format state of the original std::cout.
144 Teuchos::oblackholestream oldFormatState;
145 oldFormatState.copyfmt(std::cout);
146
147 *outStream \
148 << "===============================================================================\n" \
149 << "| |\n" \
150 << "| Unit Test (Basis_HGRAD_LINE_Cn_FEM_JACOBI) |\n" \
151 << "| |\n" \
152 << "| 1) Patch test involving mass and stiffness matrices, |\n" \
153 << "| for the Neumann problem on a REFERENCE line: |\n" \
154 << "| |\n" \
155 << "| - u'' + u = f in (-1,1), u' = g at -1,1 |\n" \
156 << "| |\n" \
157 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
158 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
159 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
160 << "| |\n" \
161 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
162 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
163 << "| |\n" \
164 << "===============================================================================\n"\
165 << "| TEST 1: Patch test |\n"\
166 << "===============================================================================\n";
167
168
169 int errorFlag = 0;
170 double zero = 100*INTREPID_TOL;
171 outStream -> precision(20);
172
173
174 try {
175
176 int max_order = 10; // max total order of polynomial solution
177
178 // Define array containing points at which the solution is evaluated
179 int numIntervals = 100;
180 int numInterpPoints = numIntervals + 1;
181 FieldContainer<double> interp_points(numInterpPoints, 1);
182 for (int i=0; i<numInterpPoints; i++) {
183 interp_points(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
184 }
185
186 DefaultCubatureFactory<double> cubFactory; // create factory
187 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
188
189 for (int soln_order=1; soln_order <= max_order; soln_order++) {
190
191 // evaluate exact solution
192 FieldContainer<double> exact_solution(1, numInterpPoints);
193 u_exact(exact_solution, interp_points, soln_order);
194
195 for (int basis_order=soln_order; basis_order <= max_order; basis_order++) {
196
197 //create basis
198 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
199 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(basis_order) );
200 int numFields = lineBasis->getCardinality();
201
202 // create cubature
203 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, 2*basis_order-2);
204 int numCubPoints = lineCub->getNumPoints();
205 int spaceDim = lineCub->getDimension();
206
207 /* Computational arrays. */
208 FieldContainer<double> cub_points(numCubPoints, spaceDim);
209 FieldContainer<double> cub_points_physical(1, numCubPoints, spaceDim);
210 FieldContainer<double> cub_weights(numCubPoints);
211 FieldContainer<double> cell_nodes(1, 2, spaceDim);
212 FieldContainer<double> jacobian(1, numCubPoints, spaceDim, spaceDim);
213 FieldContainer<double> jacobian_inv(1, numCubPoints, spaceDim, spaceDim);
214 FieldContainer<double> jacobian_det(1, numCubPoints);
215 FieldContainer<double> weighted_measure(1, numCubPoints);
216
217 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints);
218 FieldContainer<double> transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
219 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
220 FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
221 FieldContainer<double> transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
222 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
223 FieldContainer<double> fe_matrix(1, numFields, numFields);
224
225 FieldContainer<double> rhs_at_cub_points_physical(1, numCubPoints);
226 FieldContainer<double> rhs_and_soln_vector(1, numFields);
227
228 FieldContainer<double> one_point(1, 1);
229 FieldContainer<double> value_of_basis_at_one(numFields, 1);
230 FieldContainer<double> value_of_basis_at_minusone(numFields, 1);
231 FieldContainer<double> bc_neumann(2, numFields);
232
233 FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
234 FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
235 FieldContainer<double> interpolant(1, numInterpPoints);
236
237 FieldContainer<int> ipiv(numFields);
238
239 /******************* START COMPUTATION ***********************/
240
241 // get cubature points and weights
242 lineCub->getCubature(cub_points, cub_weights);
243
244 // fill cell vertex array
245 cell_nodes(0, 0, 0) = -1.0;
246 cell_nodes(0, 1, 0) = 1.0;
247
248 // compute geometric cell information
249 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, line);
250 CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
251 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
252
253 // compute weighted measure
254 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
255
257 // Computing mass matrices:
258 // tabulate values of basis functions at (reference) cubature points
259 lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
260
261 // transform values of basis functions
262 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
263 value_of_basis_at_cub_points);
264
265 // multiply with weighted measure
266 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
267 weighted_measure,
268 transformed_value_of_basis_at_cub_points);
269
270 // compute mass matrices
271 FunctionSpaceTools::integrate<double>(fe_matrix,
272 transformed_value_of_basis_at_cub_points,
273 weighted_transformed_value_of_basis_at_cub_points,
274 COMP_CPP);
276
278 // Computing stiffness matrices:
279 // tabulate gradients of basis functions at (reference) cubature points
280 lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
281
282 // transform gradients of basis functions
283 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
284 jacobian_inv,
285 grad_of_basis_at_cub_points);
286
287 // multiply with weighted measure
288 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
289 weighted_measure,
290 transformed_grad_of_basis_at_cub_points);
291
292 // compute stiffness matrices and sum into fe_matrix
293 FunctionSpaceTools::integrate<double>(fe_matrix,
294 transformed_grad_of_basis_at_cub_points,
295 weighted_transformed_grad_of_basis_at_cub_points,
296 COMP_CPP,
297 true);
299
301 // Computing RHS contributions:
302 // map (reference) cubature points to physical space
303 CellTools<double>::mapToPhysicalFrame(cub_points_physical, cub_points, cell_nodes, line);
304
305 // evaluate rhs function
306 rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
307
308 // compute rhs
309 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
310 rhs_at_cub_points_physical,
311 weighted_transformed_value_of_basis_at_cub_points,
312 COMP_CPP);
313
314 // compute neumann b.c. contributions and adjust rhs
315 one_point(0,0) = 1.0; lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
316 one_point(0,0) = -1.0; lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
317 neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
318 for (int i=0; i<numFields; i++) {
319 rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
320 rhs_and_soln_vector(0, i) += bc_neumann(1, i);
321 }
323
325 // Solution of linear system:
326 int info = 0;
327 Teuchos::LAPACK<int, double> solver;
328 //solver.GESV(numRows, 1, &fe_mat(0,0), numRows, &ipiv(0), &fe_vec(0), numRows, &info);
329 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
331
333 // Building interpolant:
334 // evaluate basis at interpolation points
335 lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
336 // transform values of basis functions
337 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
338 value_of_basis_at_interp_points);
339 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
341
342 /******************* END COMPUTATION ***********************/
343
344 RealSpaceTools<double>::subtract(interpolant, exact_solution);
345
346 *outStream << "\nNorm-2 difference between exact solution polynomial of order "
347 << soln_order << " and finite element interpolant of order " << basis_order << ": "
348 << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) << "\n";
349
350 if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) > zero) {
351 *outStream << "\n\nPatch test failed for solution polynomial order "
352 << soln_order << " and basis order " << basis_order << "\n\n";
353 errorFlag++;
354 }
355
356 } // end for basis_order
357
358 } // end for soln_order
359
360 }
361 // Catch unexpected errors
362 catch (const std::logic_error & err) {
363 *outStream << err.what() << "\n\n";
364 errorFlag = -1000;
365 };
366
367 if (errorFlag != 0)
368 std::cout << "End Result: TEST FAILED\n";
369 else
370 std::cout << "End Result: TEST PASSED\n";
371
372 // reset format state of std::cout
373 std::cout.copyfmt(oldFormatState);
374
375 return errorFlag;
376}
void u_exact(FieldContainer< double > &, const FieldContainer< double > &, int)
exact solution
Definition: test_02.cpp:118
void neumann(FieldContainer< double > &, const FieldContainer< double > &, const FieldContainer< double > &, int)
neumann boundary conditions
Definition: test_02.cpp:96
void rhsFunc(FieldContainer< double > &, const FieldContainer< double > &, int)
right-hand side function
Definition: test_02.cpp:71
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_LINE_Cn_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
A stateless class for operations on cell data. Provides methods for:
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int dimension(const int whichDim) const
Returns the specified dimension.
Implementation of basic linear algebra functionality in Euclidean space.
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2.
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.