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 
53 #include "Intrepid_ArrayTools.hpp"
55 #include "Intrepid_CellTools.hpp"
56 #include "Teuchos_oblackholestream.hpp"
57 #include "Teuchos_RCP.hpp"
58 #include "Teuchos_GlobalMPISession.hpp"
59 #include "Teuchos_SerialDenseMatrix.hpp"
60 #include "Teuchos_SerialDenseVector.hpp"
61 #include "Teuchos_LAPACK.hpp"
62 
63 using namespace std;
64 using namespace Intrepid;
65 
66 void rhsFunc(FieldContainer<double> &, const FieldContainer<double> &, int, int);
67 void neumann(FieldContainer<double> & ,
68  const FieldContainer<double> & ,
69  const FieldContainer<double> & ,
70  const shards::CellTopology & ,
71  int, int, int);
72 void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int, int);
73 
75 void rhsFunc(FieldContainer<double> & result,
76  const FieldContainer<double> & points,
77  int xd,
78  int yd) {
79 
80  int x = 0, y = 1;
81 
82  // second x-derivatives of u
83  if (xd > 1) {
84  for (int cell=0; cell<result.dimension(0); cell++) {
85  for (int pt=0; pt<result.dimension(1); pt++) {
86  result(cell,pt) = - xd*(xd-1)*std::pow(points(cell,pt,x), xd-2) * std::pow(points(cell,pt,y), yd);
87  }
88  }
89  }
90 
91  // second y-derivatives of u
92  if (yd > 1) {
93  for (int cell=0; cell<result.dimension(0); cell++) {
94  for (int pt=0; pt<result.dimension(1); pt++) {
95  result(cell,pt) -= yd*(yd-1)*std::pow(points(cell,pt,y), yd-2) * std::pow(points(cell,pt,x), xd);
96  }
97  }
98  }
99 
100  // add u
101  for (int cell=0; cell<result.dimension(0); cell++) {
102  for (int pt=0; pt<result.dimension(1); pt++) {
103  result(cell,pt) += std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
104  }
105  }
106 
107 }
108 
109 
111 void neumann(FieldContainer<double> & result,
112  const FieldContainer<double> & points,
113  const FieldContainer<double> & jacs,
114  const shards::CellTopology & parentCell,
115  int sideOrdinal, int xd, int yd) {
116 
117  int x = 0, y = 1;
118 
119  int numCells = result.dimension(0);
120  int numPoints = result.dimension(1);
121 
122  FieldContainer<double> grad_u(numCells, numPoints, 2);
123  FieldContainer<double> side_normals(numCells, numPoints, 2);
124  FieldContainer<double> normal_lengths(numCells, numPoints);
125 
126  // first x-derivatives of u
127  if (xd > 0) {
128  for (int cell=0; cell<numCells; cell++) {
129  for (int pt=0; pt<numPoints; pt++) {
130  grad_u(cell,pt,x) = xd*std::pow(points(cell,pt,x), xd-1) * std::pow(points(cell,pt,y), yd);
131  }
132  }
133  }
134 
135  // first y-derivatives of u
136  if (yd > 0) {
137  for (int cell=0; cell<numCells; cell++) {
138  for (int pt=0; pt<numPoints; pt++) {
139  grad_u(cell,pt,y) = yd*std::pow(points(cell,pt,y), yd-1) * std::pow(points(cell,pt,x), xd);
140  }
141  }
142  }
143 
144  CellTools<double>::getPhysicalSideNormals(side_normals, jacs, sideOrdinal, parentCell);
145 
146  // scale normals
147  RealSpaceTools<double>::vectorNorm(normal_lengths, side_normals, NORM_TWO);
148  FunctionSpaceTools::scalarMultiplyDataData<double>(side_normals, normal_lengths, side_normals, true);
149 
150  FunctionSpaceTools::dotMultiplyDataData<double>(result, grad_u, side_normals);
151 
152 }
153 
155 void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int xd, int yd) {
156  int x = 0, y = 1;
157  for (int cell=0; cell<result.dimension(0); cell++) {
158  for (int pt=0; pt<result.dimension(1); pt++) {
159  result(cell,pt) = std::pow(points(pt,x), xd)*std::pow(points(pt,y), yd);
160  }
161  }
162 }
163 
164 
165 
166 
167 int main(int argc, char *argv[]) {
168 
169  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
170 
171  // This little trick lets us print to std::cout only if
172  // a (dummy) command-line argument is provided.
173  int iprint = argc - 1;
174  Teuchos::RCP<std::ostream> outStream;
175  Teuchos::oblackholestream bhs; // outputs nothing
176  if (iprint > 0)
177  outStream = Teuchos::rcp(&std::cout, false);
178  else
179  outStream = Teuchos::rcp(&bhs, false);
180 
181  // Save the format state of the original std::cout.
182  Teuchos::oblackholestream oldFormatState;
183  oldFormatState.copyfmt(std::cout);
184 
185  *outStream \
186  << "===============================================================================\n" \
187  << "| |\n" \
188  << "| Unit Test (Basis_HGRAD_TRI_C2_FEM) |\n" \
189  << "| |\n" \
190  << "| 1) Patch test involving mass and stiffness matrices, |\n" \
191  << "| for the Neumann problem on a triangular patch |\n" \
192  << "| Omega with boundary Gamma. |\n" \
193  << "| |\n" \
194  << "| - div (grad u) + u = f in Omega, (grad u) . n = g on Gamma |\n" \
195  << "| |\n" \
196  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
197  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
198  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
199  << "| |\n" \
200  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
201  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
202  << "| |\n" \
203  << "===============================================================================\n"\
204  << "| TEST 1: Patch test |\n"\
205  << "===============================================================================\n";
206 
207 
208  int errorFlag = 0;
209 
210  outStream -> precision(16);
211 
212 
213  try {
214 
215  int max_order = 2; // max total order of polynomial solution
216  DefaultCubatureFactory<double> cubFactory; // create cubature factory
217  shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >()); // create parent cell topology
218  shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >()); // create relevant subcell (side) topology
219  int cellDim = cell.getDimension();
220  int sideDim = side.getDimension();
221 
222  // Define array containing points at which the solution is evaluated, in reference cell.
223  int numIntervals = 10;
224  int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2;
225  FieldContainer<double> interp_points_ref(numInterpPoints, 2);
226  int counter = 0;
227  for (int j=0; j<=numIntervals; j++) {
228  for (int i=0; i<=numIntervals; i++) {
229  if (i <= numIntervals-j) {
230  interp_points_ref(counter,0) = i*(1.0/numIntervals);
231  interp_points_ref(counter,1) = j*(1.0/numIntervals);
232  counter++;
233  }
234  }
235  }
236 
237  /* Parent cell definition. */
238  FieldContainer<double> cell_nodes(1, 3, cellDim);
239  // Perturbed reference triangle.
240  cell_nodes(0, 0, 0) = 0.1;
241  cell_nodes(0, 0, 1) = -0.1;
242  cell_nodes(0, 1, 0) = 1.1;
243  cell_nodes(0, 1, 1) = -0.1;
244  cell_nodes(0, 2, 0) = 0.1;
245  cell_nodes(0, 2, 1) = 0.9;
246  // Reference triangle.
247  /*cell_nodes(0, 0, 0) = 0.0;
248  cell_nodes(0, 0, 1) = 0.0;
249  cell_nodes(0, 1, 0) = 1.0;
250  cell_nodes(0, 1, 1) = 0.0;
251  cell_nodes(0, 2, 0) = 0.0;
252  cell_nodes(0, 2, 1) = 1.0;*/
253 
254  FieldContainer<double> interp_points(1, numInterpPoints, cellDim);
255  CellTools<double>::mapToPhysicalFrame(interp_points, interp_points_ref, cell_nodes, cell);
256  interp_points.resize(numInterpPoints, cellDim);
257 
258  for (int x_order=0; x_order <= max_order; x_order++) {
259  for (int y_order=0; y_order <= max_order-x_order; y_order++) {
260 
261  // evaluate exact solution
262  FieldContainer<double> exact_solution(1, numInterpPoints);
263  u_exact(exact_solution, interp_points, x_order, y_order);
264 
265  int basis_order = 2;
266 
267  // set test tolerance
268  double zero = basis_order*basis_order*100*INTREPID_TOL;
269 
270  //create basis
271  Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
272  Teuchos::rcp(new Basis_HGRAD_TRI_C2_FEM<double,FieldContainer<double> >() );
273  int numFields = basis->getCardinality();
274 
275  // create cubatures
276  Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
277  Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*basis_order);
278  int numCubPointsCell = cellCub->getNumPoints();
279  int numCubPointsSide = sideCub->getNumPoints();
280 
281  /* Computational arrays. */
282  /* Section 1: Related to parent cell integration. */
283  FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
284  FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
285  FieldContainer<double> cub_weights_cell(numCubPointsCell);
286  FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
287  FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
288  FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
289  FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
290 
291  FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
292  FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
293  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
294  FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
295  FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
296  FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
297  FieldContainer<double> fe_matrix(1, numFields, numFields);
298 
299  FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
300  FieldContainer<double> rhs_and_soln_vector(1, numFields);
301 
302  /* Section 2: Related to subcell (side) integration. */
303  unsigned numSides = 3;
304  FieldContainer<double> cub_points_side(numCubPointsSide, sideDim);
305  FieldContainer<double> cub_weights_side(numCubPointsSide);
306  FieldContainer<double> cub_points_side_refcell(numCubPointsSide, cellDim);
307  FieldContainer<double> cub_points_side_physical(1, numCubPointsSide, cellDim);
308  FieldContainer<double> jacobian_side_refcell(1, numCubPointsSide, cellDim, cellDim);
309  FieldContainer<double> jacobian_det_side_refcell(1, numCubPointsSide);
310  FieldContainer<double> weighted_measure_side_refcell(1, numCubPointsSide);
311 
312  FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields, numCubPointsSide);
313  FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
314  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
315  FieldContainer<double> neumann_data_at_cub_points_side_physical(1, numCubPointsSide);
316  FieldContainer<double> neumann_fields_per_side(1, numFields);
317 
318  /* Section 3: Related to global interpolant. */
319  FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
320  FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
321  FieldContainer<double> interpolant(1, numInterpPoints);
322 
323  FieldContainer<int> ipiv(numFields);
324 
325 
326 
327  /******************* START COMPUTATION ***********************/
328 
329  // get cubature points and weights
330  cellCub->getCubature(cub_points_cell, cub_weights_cell);
331 
332  // compute geometric cell information
333  CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes, cell);
334  CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
335  CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
336 
337  // compute weighted measure
338  FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
339 
341  // Computing mass matrices:
342  // tabulate values of basis functions at (reference) cubature points
343  basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
344 
345  // transform values of basis functions
346  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
347  value_of_basis_at_cub_points_cell);
348 
349  // multiply with weighted measure
350  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
351  weighted_measure_cell,
352  transformed_value_of_basis_at_cub_points_cell);
353 
354  // compute mass matrices
355  FunctionSpaceTools::integrate<double>(fe_matrix,
356  transformed_value_of_basis_at_cub_points_cell,
357  weighted_transformed_value_of_basis_at_cub_points_cell,
358  COMP_BLAS);
360 
362  // Computing stiffness matrices:
363  // tabulate gradients of basis functions at (reference) cubature points
364  basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
365 
366  // transform gradients of basis functions
367  FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
368  jacobian_inv_cell,
369  grad_of_basis_at_cub_points_cell);
370 
371  // multiply with weighted measure
372  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
373  weighted_measure_cell,
374  transformed_grad_of_basis_at_cub_points_cell);
375 
376  // compute stiffness matrices and sum into fe_matrix
377  FunctionSpaceTools::integrate<double>(fe_matrix,
378  transformed_grad_of_basis_at_cub_points_cell,
379  weighted_transformed_grad_of_basis_at_cub_points_cell,
380  COMP_BLAS,
381  true);
383 
385  // Computing RHS contributions:
386  // map cell (reference) cubature points to physical space
387  CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes, cell);
388 
389  // evaluate rhs function
390  rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order);
391 
392  // compute rhs
393  FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
394  rhs_at_cub_points_cell_physical,
395  weighted_transformed_value_of_basis_at_cub_points_cell,
396  COMP_BLAS);
397 
398  // compute neumann b.c. contributions and adjust rhs
399  sideCub->getCubature(cub_points_side, cub_weights_side);
400  for (unsigned i=0; i<numSides; i++) {
401  // compute geometric cell information
402  CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
403  CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes, cell);
404  CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
405 
406  // compute weighted edge measure
407  FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
408  jacobian_side_refcell,
409  cub_weights_side,
410  i,
411  cell);
412 
413  // tabulate values of basis functions at side cubature points, in the reference parent cell domain
414  basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
415  // transform
416  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
417  value_of_basis_at_cub_points_side_refcell);
418 
419  // multiply with weighted measure
420  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
421  weighted_measure_side_refcell,
422  transformed_value_of_basis_at_cub_points_side_refcell);
423 
424  // compute Neumann data
425  // map side cubature points in reference parent cell domain to physical space
426  CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes, cell);
427  // now compute data
428  neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
429  cell, (int)i, x_order, y_order);
430 
431  FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
432  neumann_data_at_cub_points_side_physical,
433  weighted_transformed_value_of_basis_at_cub_points_side_refcell,
434  COMP_BLAS);
435 
436  // adjust RHS
437  RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
438  }
440 
442  // Solution of linear system:
443  int info = 0;
444  Teuchos::LAPACK<int, double> solver;
445  solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
447 
449  // Building interpolant:
450  // evaluate basis at interpolation points
451  basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
452  // transform values of basis functions
453  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
454  value_of_basis_at_interp_points);
455  FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
457 
458  /******************* END COMPUTATION ***********************/
459 
460  RealSpaceTools<double>::subtract(interpolant, exact_solution);
461 
462  *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
463  << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
464  << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
465  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
466 
467  if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
468  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
469  *outStream << "\n\nPatch test failed for solution polynomial order ("
470  << x_order << ", " << y_order << ") and basis order " << basis_order << "\n\n";
471  errorFlag++;
472  }
473  } // end for y_order
474  } // end for x_order
475 
476  }
477  // Catch unexpected errors
478  catch (const std::logic_error & err) {
479  *outStream << err.what() << "\n\n";
480  errorFlag = -1000;
481  };
482 
483  if (errorFlag != 0)
484  std::cout << "End Result: TEST FAILED\n";
485  else
486  std::cout << "End Result: TEST PASSED\n";
487 
488  // reset format state of std::cout
489  std::cout.copyfmt(oldFormatState);
490 
491  return errorFlag;
492 }
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 Intrepid::HGRAD_TRI_C2_FEM class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
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.
A stateless class for operations on cell data. Provides methods for:
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell...