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_C1_FEM) |\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 = 1; // 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 //create basis
190 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
191 Teuchos::rcp(new Basis_HGRAD_LINE_C1_FEM<double,FieldContainer<double> >() );
192 int numFields = lineBasis->getCardinality();
193 int basis_order = lineBasis->getDegree();
194
195 // create cubature
196 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, 2);
197 int numCubPoints = lineCub->getNumPoints();
198 int spaceDim = lineCub->getDimension();
199
200 for (int soln_order=0; soln_order <= max_order; soln_order++) {
201
202 // evaluate exact solution
203 FieldContainer<double> exact_solution(1, numInterpPoints);
204 u_exact(exact_solution, interp_points, soln_order);
205
206 /* Computational arrays. */
207 FieldContainer<double> cub_points(numCubPoints, spaceDim);
208 FieldContainer<double> cub_points_physical(1, numCubPoints, spaceDim);
209 FieldContainer<double> cub_weights(numCubPoints);
210 FieldContainer<double> cell_nodes(1, 2, spaceDim);
211 FieldContainer<double> jacobian(1, numCubPoints, spaceDim, spaceDim);
212 FieldContainer<double> jacobian_inv(1, numCubPoints, spaceDim, spaceDim);
213 FieldContainer<double> jacobian_det(1, numCubPoints);
214 FieldContainer<double> weighted_measure(1, numCubPoints);
215
216 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints);
217 FieldContainer<double> transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
218 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
219 FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
220 FieldContainer<double> transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
221 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
222 FieldContainer<double> fe_matrix(1, numFields, numFields);
223
224 FieldContainer<double> rhs_at_cub_points_physical(1, numCubPoints);
225 FieldContainer<double> rhs_and_soln_vector(1, numFields);
226
227 FieldContainer<double> one_point(1, 1);
228 FieldContainer<double> value_of_basis_at_one(numFields, 1);
229 FieldContainer<double> value_of_basis_at_minusone(numFields, 1);
230 FieldContainer<double> bc_neumann(2, numFields);
231
232 FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
233 FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
234 FieldContainer<double> interpolant(1, numInterpPoints);
235
236 FieldContainer<int> ipiv(numFields);
237
238 /******************* START COMPUTATION ***********************/
239
240 // get cubature points and weights
241 lineCub->getCubature(cub_points, cub_weights);
242
243 // fill cell vertex array
244 cell_nodes(0, 0, 0) = -1.0;
245 cell_nodes(0, 1, 0) = 1.0;
246
247 // compute geometric cell information
248 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, line);
249 CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
250 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
251
252 // compute weighted measure
253 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
254
256 // Computing mass matrices:
257 // tabulate values of basis functions at (reference) cubature points
258 lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
259
260 // transform values of basis functions
261 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
262 value_of_basis_at_cub_points);
263
264 // multiply with weighted measure
265 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
266 weighted_measure,
267 transformed_value_of_basis_at_cub_points);
268
269 // compute mass matrices
270 FunctionSpaceTools::integrate<double>(fe_matrix,
271 transformed_value_of_basis_at_cub_points,
272 weighted_transformed_value_of_basis_at_cub_points,
273 COMP_CPP);
275
277 // Computing stiffness matrices:
278 // tabulate gradients of basis functions at (reference) cubature points
279 lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
280
281 // transform gradients of basis functions
282 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
283 jacobian_inv,
284 grad_of_basis_at_cub_points);
285
286 // multiply with weighted measure
287 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
288 weighted_measure,
289 transformed_grad_of_basis_at_cub_points);
290
291 // compute stiffness matrices and sum into fe_matrix
292 FunctionSpaceTools::integrate<double>(fe_matrix,
293 transformed_grad_of_basis_at_cub_points,
294 weighted_transformed_grad_of_basis_at_cub_points,
295 COMP_CPP,
296 true);
298
300 // Computing RHS contributions:
301 // map (reference) cubature points to physical space
302 CellTools<double>::mapToPhysicalFrame(cub_points_physical, cub_points, cell_nodes, line);
303
304 // evaluate rhs function
305 rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
306
307 // compute rhs
308 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
309 rhs_at_cub_points_physical,
310 weighted_transformed_value_of_basis_at_cub_points,
311 COMP_CPP);
312
313 // compute neumann b.c. contributions and adjust rhs
314 one_point(0,0) = 1.0; lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
315 one_point(0,0) = -1.0; lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
316 neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
317 for (int i=0; i<numFields; i++) {
318 rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
319 rhs_and_soln_vector(0, i) += bc_neumann(1, i);
320 }
322
324 // Solution of linear system:
325 int info = 0;
326 Teuchos::LAPACK<int, double> solver;
327 //solver.GESV(numRows, 1, &fe_mat(0,0), numRows, &ipiv(0), &fe_vec(0), numRows, &info);
328 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
330
332 // Building interpolant:
333 // evaluate basis at interpolation points
334 lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
335 // transform values of basis functions
336 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
337 value_of_basis_at_interp_points);
338 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
340
341 /******************* END COMPUTATION ***********************/
342
343 RealSpaceTools<double>::subtract(interpolant, exact_solution);
344
345 *outStream << "\nNorm-2 difference between exact solution polynomial of order "
346 << soln_order << " and finite element interpolant of order " << basis_order << ": "
347 << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) << "\n";
348
349 if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) > zero) {
350 *outStream << "\n\nPatch test failed for solution polynomial order "
351 << soln_order << " and basis order " << basis_order << "\n\n";
352 errorFlag++;
353 }
354
355 } // end for soln_order
356
357 }
358 // Catch unexpected errors
359 catch (const std::logic_error & err) {
360 *outStream << err.what() << "\n\n";
361 errorFlag = -1000;
362 };
363
364 if (errorFlag != 0)
365 std::cout << "End Result: TEST FAILED\n";
366 else
367 std::cout << "End Result: TEST PASSED\n";
368
369 // reset format state of std::cout
370 std::cout.copyfmt(oldFormatState);
371
372 return errorFlag;
373}
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_C1_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Line cell.
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.