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