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