Intrepid
test_02.cpp
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_HEX_Cn_FEM.hpp"
51 #include "Intrepid_HDIV_HEX_In_FEM.hpp"
52 #include "Intrepid_DefaultCubatureFactory.hpp"
53 #include "Intrepid_PointTools.hpp"
54 #include "Intrepid_RealSpaceTools.hpp"
55 #include "Intrepid_ArrayTools.hpp"
56 #include "Intrepid_FunctionSpaceTools.hpp"
57 #include "Intrepid_CellTools.hpp"
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 
65 using namespace std;
66 using namespace Intrepid;
67 
68 void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int );
69 void 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
73 void rhsFunc( FieldContainer<double> &result,
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 
99 void u_exact( FieldContainer<double> &result,
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 
114 int 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
365  CellTools<double>::getReferenceSideNormal(side_normal ,
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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Intrepid::Basis_HDIV_HEX_In_FEM
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Definition: Intrepid_HDIV_HEX_In_FEM.hpp:70
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Intrepid_HDIV_HEX_In_FEM.hpp
Header file for the Intrepid::HDIV_HEX_In_FEM class.
Intrepid_PointTools.hpp
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
Intrepid::FieldContainer::dimension
int dimension(const int whichDim) const
Returns the specified dimension.
Definition: Intrepid_FieldContainerDef.hpp:658
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid_HGRAD_HEX_Cn_FEM.hpp
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::Basis_HGRAD_HEX_Cn_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Definition: Intrepid_HGRAD_HEX_Cn_FEM.hpp:72
Intrepid_FunctionSpaceTools.hpp
Header file for the Intrepid::FunctionSpaceTools class.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
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::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:109
Generated by   doxygen 1.8.5