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
58#include "Teuchos_oblackholestream.hpp"
59#include "Teuchos_RCP.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
61#include "Teuchos_SerialDenseMatrix.hpp"
62#include "Teuchos_SerialDenseVector.hpp"
63#include "Teuchos_LAPACK.hpp"
64
65using namespace std;
66using namespace Intrepid;
67
68void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int );
69void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int, int, int );
70
71// This is the rhs for (div tau,w) = (f,w),
72// which makes f the negative Laplacian of scalar solution
74 const FieldContainer<double> &points,
75 int xd,
76 int yd ,
77 int zd)
78{
79 for (int cell=0;cell<result.dimension(0);cell++) {
80 for (int pt=0;pt<result.dimension(1);pt++) {
81 result(cell,pt) = 0.0;
82 if (xd >=2) {
83 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd)
84 *pow(points(cell,pt,2),zd);
85 }
86 if (yd >=2) {
87 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2)
88 *pow(points(cell,pt,2),zd);
89 }
90 if (zd>=2) {
91 result(cell,pt) += zd*(zd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd)
92 *pow(points(cell,pt,2),zd-2);
93
94 }
95 }
96 }
97}
98
100 const FieldContainer<double> &points,
101 int xd,
102 int yd,
103 int zd)
104{
105 for (int cell=0;cell<result.dimension(0);cell++){
106 for (int pt=0;pt<result.dimension(1);pt++) {
107 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
108 *std::pow(points(cell,pt,2),zd);
109 }
110 }
111 return;
112}
113
114int main(int argc, char *argv[]) {
115 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
116
117 // This little trick lets us print to std::cout only if
118 // a (dummy) command-line argument is provided.
119 int iprint = argc - 1;
120 Teuchos::RCP<std::ostream> outStream;
121 Teuchos::oblackholestream bhs; // outputs nothing
122 if (iprint > 0)
123 outStream = Teuchos::rcp(&std::cout, false);
124 else
125 outStream = Teuchos::rcp(&bhs, false);
126
127 // Save the format state of the original std::cout.
128 Teuchos::oblackholestream oldFormatState;
129 oldFormatState.copyfmt(std::cout);
130
131 *outStream \
132 << "===============================================================================\n" \
133 << "| |\n" \
134 << "| Unit Test (Basis_HGRAD_HEX_In_FEM) |\n" \
135 << "| |\n" \
136 << "| 1) Patch test involving H(div) matrices |\n" \
137 << "| for the Dirichlet problem on a hexahedron |\n" \
138 << "| Omega with boundary Gamma. |\n" \
139 << "| |\n" \
140 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
141 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
142 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
143 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
144 << "| |\n" \
145 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
146 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
147 << "| |\n" \
148 << "===============================================================================\n" \
149 << "| TEST 1: Patch test |\n" \
150 << "===============================================================================\n";
151
152
153 int errorFlag = 0;
154
155 outStream -> precision(16);
156
157 try {
158 DefaultCubatureFactory<double> cubFactory; // create cubature factory
159 shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<> >()); // create parent cell topology
160 shards::CellTopology side(shards::getCellTopologyData< shards::Quadrilateral<> >()); // create relevant subcell (side) topology
161 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >() ); // for getting points to construct the basis
162
163 int cellDim = cell.getDimension();
164 int sideDim = side.getDimension();
165
166 int min_order = 0;
167 int max_order = 3;
168
169 int numIntervals = 2;
170 int numInterpPoints = (numIntervals + 1)*(numIntervals + 1)*(numIntervals+1);
171 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
172 int counter = 0;
173 for (int k=0;k<numIntervals;k++) {
174 for (int j=0; j<=numIntervals; j++) {
175 for (int i=0; i<=numIntervals; i++) {
176 interp_points_ref(counter,0) = i*(1.0/numIntervals);
177 interp_points_ref(counter,1) = j*(1.0/numIntervals);
178 interp_points_ref(counter,2) = k*(1.0/numIntervals);
179 counter++;
180 }
181 }
182 }
183
184 for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
185 // create bases
186 // get the points for the vector basis
187 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
188 Teuchos::rcp(new Basis_HDIV_HEX_In_FEM<double,FieldContainer<double> >(basis_order+1,POINTTYPE_SPECTRAL) );
189
190 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
191 Teuchos::rcp(new Basis_HGRAD_HEX_Cn_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_SPECTRAL) );
192
193 int numVectorFields = vectorBasis->getCardinality();
194 int numScalarFields = scalarBasis->getCardinality();
195 int numTotalFields = numVectorFields + numScalarFields;
196
197 // create cubatures
198 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
199 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
200
201 int numCubPointsCell = cellCub->getNumPoints();
202 int numCubPointsSide = sideCub->getNumPoints();
203
204 // hold cubature information
205 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
206 FieldContainer<double> cub_weights_cell(numCubPointsCell);
207 FieldContainer<double> cub_points_side( numCubPointsSide, sideDim );
208 FieldContainer<double> cub_weights_side( numCubPointsSide );
209 FieldContainer<double> cub_points_side_refcell( numCubPointsSide , cellDim );
210
211 // hold basis function information on refcell
212 FieldContainer<double> value_of_v_basis_at_cub_points_cell(numVectorFields, numCubPointsCell, cellDim );
213 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
214 FieldContainer<double> div_of_v_basis_at_cub_points_cell( numVectorFields, numCubPointsCell );
215 FieldContainer<double> w_div_of_v_basis_at_cub_points_cell( 1, numVectorFields , numCubPointsCell );
216 FieldContainer<double> value_of_s_basis_at_cub_points_cell(numScalarFields,numCubPointsCell);
217 FieldContainer<double> w_value_of_s_basis_at_cub_points_cell(1,numScalarFields,numCubPointsCell);
218
219 // containers for side integration:
220 // I just need the normal component of the vector basis
221 // and the exact solution at the cub points
222 FieldContainer<double> value_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide,cellDim);
223 FieldContainer<double> n_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide);
224 FieldContainer<double> w_n_of_v_basis_at_cub_points_side(1,numVectorFields,numCubPointsSide);
225 FieldContainer<double> diri_data_at_cub_points_side(1,numCubPointsSide);
226 FieldContainer<double> side_normal(cellDim);
227
228 // holds rhs data
229 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell);
230
231 // FEM matrices and vectors
232 FieldContainer<double> fe_matrix_M(1,numVectorFields,numVectorFields);
233 FieldContainer<double> fe_matrix_B(1,numVectorFields,numScalarFields);
234 FieldContainer<double> fe_matrix(1,numTotalFields,numTotalFields);
235
236 FieldContainer<double> rhs_vector_vec(1,numVectorFields);
237 FieldContainer<double> rhs_vector_scal(1,numScalarFields);
238 FieldContainer<double> rhs_and_soln_vec(1,numTotalFields);
239
240 FieldContainer<int> ipiv(numTotalFields);
241 FieldContainer<double> value_of_s_basis_at_interp_points( numScalarFields , numInterpPoints);
242 FieldContainer<double> interpolant( 1 , numInterpPoints );
243
244 // set test tolerance
245 double zero = (basis_order+1)*(basis_order+1)*1000.0*INTREPID_TOL;
246
247 // build matrices outside the loop, and then just do the rhs
248 // for each iteration
249
250 cellCub->getCubature(cub_points_cell, cub_weights_cell);
251 sideCub->getCubature(cub_points_side, cub_weights_side);
252
253 // need the vector basis & its divergences
254 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
255 cub_points_cell,
256 OPERATOR_VALUE);
257 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
258 cub_points_cell,
259 OPERATOR_DIV);
260
261 // need the scalar basis as well
262 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
263 cub_points_cell,
264 OPERATOR_VALUE);
265
266 // construct mass matrix
267 cub_weights_cell.resize(1,numCubPointsCell);
268 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
269 cub_weights_cell ,
270 value_of_v_basis_at_cub_points_cell );
271 cub_weights_cell.resize(numCubPointsCell);
272
273
274 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
275 FunctionSpaceTools::integrate<double>(fe_matrix_M,
276 w_value_of_v_basis_at_cub_points_cell ,
277 value_of_v_basis_at_cub_points_cell ,
278 COMP_BLAS );
279 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
280
281 // div matrix
282 cub_weights_cell.resize(1,numCubPointsCell);
283 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
284 cub_weights_cell,
285 div_of_v_basis_at_cub_points_cell);
286 cub_weights_cell.resize(numCubPointsCell);
287
288 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
289 FunctionSpaceTools::integrate<double>(fe_matrix_B,
290 w_div_of_v_basis_at_cub_points_cell ,
291 value_of_s_basis_at_cub_points_cell ,
292 COMP_BLAS );
293 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
294
295 for (int x_order=0;x_order<=basis_order;x_order++) {
296 for (int y_order=0;y_order<=basis_order;y_order++) {
297 for (int z_order=0;z_order<=basis_order;z_order++) {
298
299
300 // reset global matrix since I destroyed it in LU factorization.
301 fe_matrix.initialize();
302 // insert mass matrix into global matrix
303 for (int i=0;i<numVectorFields;i++) {
304 for (int j=0;j<numVectorFields;j++) {
305 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
306 }
307 }
308
309 // insert div matrix into global matrix
310 for (int i=0;i<numVectorFields;i++) {
311 for (int j=0;j<numScalarFields;j++) {
312 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
313 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
314 }
315 }
316
317 // clear old vector data
318 rhs_vector_vec.initialize();
319 rhs_vector_scal.initialize();
320 rhs_and_soln_vec.initialize();
321
322 // now get rhs vector
323 // rhs_vector_scal is just (rhs,w) for w in the scalar basis
324 // I already have the scalar basis tabulated.
325 cub_points_cell.resize(1,numCubPointsCell,cellDim);
326 rhsFunc(rhs_at_cub_points_cell,
327 cub_points_cell,
328 x_order,
329 y_order,
330 z_order);
331
332 cub_points_cell.resize(numCubPointsCell,cellDim);
333
334 cub_weights_cell.resize(1,numCubPointsCell);
335 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
336 cub_weights_cell,
337 value_of_s_basis_at_cub_points_cell);
338 cub_weights_cell.resize(numCubPointsCell);
339 FunctionSpaceTools::integrate<double>(rhs_vector_scal,
340 rhs_at_cub_points_cell,
341 w_value_of_s_basis_at_cub_points_cell,
342 COMP_BLAS);
343
344 for (int i=0;i<numScalarFields;i++) {
345 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
346 }
347
348
349 // now get <u,v.n> on boundary
350 for (unsigned side_cur=0;side_cur<6;side_cur++) {
351 // map side cubature to current side
352 CellTools<double>::mapToReferenceSubcell( cub_points_side_refcell ,
353 cub_points_side ,
354 sideDim ,
355 (int)side_cur ,
356 cell );
357 // Evaluate dirichlet data
358 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
359 u_exact(diri_data_at_cub_points_side,
360 cub_points_side_refcell,x_order,y_order,z_order);
361
362 cub_points_side_refcell.resize(numCubPointsSide,cellDim);
363
364 // get normal direction, this has the edge weight factored into it already
366 (int)side_cur,cell );
367
368 // v.n at cub points on side
369 vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
370 cub_points_side_refcell ,
371 OPERATOR_VALUE );
372
373 for (int i=0;i<numVectorFields;i++) {
374 for (int j=0;j<numCubPointsSide;j++) {
375 n_of_v_basis_at_cub_points_side(i,j) = 0.0;
376 for (int k=0;k<cellDim;k++) {
377 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) *
378 value_of_v_basis_at_cub_points_side(i,j,k);
379 }
380 }
381 }
382
383 cub_weights_side.resize(1,numCubPointsSide);
384 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
385 cub_weights_side,
386 n_of_v_basis_at_cub_points_side);
387 cub_weights_side.resize(numCubPointsSide);
388
389 FunctionSpaceTools::integrate<double>(rhs_vector_vec,
390 diri_data_at_cub_points_side,
391 w_n_of_v_basis_at_cub_points_side,
392 COMP_BLAS,
393 false);
394
395 for (int i=0;i<numVectorFields;i++) {
396 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
397 }
398
399 }
400
401 // solve linear system
402 int info = 0;
403 Teuchos::LAPACK<int, double> solver;
404 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0],
405 numTotalFields, &info);
406
407 // compute interpolant; the scalar entries are last
408 scalarBasis->getValues(value_of_s_basis_at_interp_points,
409 interp_points_ref,
410 OPERATOR_VALUE);
411 for (int pt=0;pt<numInterpPoints;pt++) {
412 interpolant(0,pt)=0.0;
413 for (int i=0;i<numScalarFields;i++) {
414 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
415 * value_of_s_basis_at_interp_points(i,pt);
416 }
417 }
418
419 interp_points_ref.resize(1,numInterpPoints,cellDim);
420 // get exact solution for comparison
421 FieldContainer<double> exact_solution(1,numInterpPoints);
422 u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
423 interp_points_ref.resize(numInterpPoints,cellDim);
424
425 RealSpaceTools<double>::add(interpolant,exact_solution);
426
427 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
428
429 *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
430 << x_order << ", " << y_order << ", " << z_order
431 << ") and finite element interpolant of order " << basis_order << ": "
432 << nrm << "\n";
433
434 if (nrm > zero) {
435 *outStream << "\n\nPatch test failed for solution polynomial order ("
436 << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector) ("
437 << basis_order << ", " << basis_order+1 << ")\n\n";
438 errorFlag++;
439 }
440
441 }
442 }
443 }
444 }
445
446 }
447
448 catch (const std::logic_error & err) {
449 *outStream << err.what() << "\n\n";
450 errorFlag = -1000;
451 };
452
453 if (errorFlag != 0)
454 std::cout << "End Result: TEST FAILED\n";
455 else
456 std::cout << "End Result: TEST PASSED\n";
457
458 // reset format state of std::cout
459 std::cout.copyfmt(oldFormatState);
460
461 return errorFlag;
462}
void rhsFunc(FieldContainer< double > &, const FieldContainer< double > &, int, int, int)
right-hand side function
Definition: test_02.cpp:73
void u_exact(FieldContainer< double > &, const FieldContainer< double > &, int, int, int)
exact solution
Definition: test_02.cpp:99
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::HDIV_HEX_In_FEM class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getReferenceSideNormal(ArraySideNormal &refSideNormal, const int &sideOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
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.
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.
static void add(Scalar *sumArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Adds contiguous data inArray1 and inArray2 of size size: sumArray = inArray1 + inArray2.