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 
50 #include "Intrepid_HGRAD_TET_Cn_FEM_ORTH.hpp"
54 #include "Intrepid_ArrayTools.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
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 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for the Intrepid::CellTools class.
int dimension(const int whichDim) const
Returns the specified dimension.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(div)-compatible Raviart-Thomas basis of arbitrary degree on Tetrahedr...
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
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 the default H(grad)-compatible orthogonal basis of arbitrary degree on tetrahedron...
Header file for the Intrepid::HDIV_TET_In_FEM class.
A stateless class for operations on cell data. Provides methods for: