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_Cn_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 = 10; // 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  // we'll test two types of bases
259  EPointType pointtype[] = {POINTTYPE_EQUISPACED, POINTTYPE_WARPBLEND};
260  for (int ptype=0; ptype < 2; ptype++) {
261 
262  *outStream << "\nTesting bases with " << EPointTypeToString(pointtype[ptype]) << ":\n";
263 
264  for (int x_order=0; x_order <= max_order; x_order++) {
265  for (int y_order=0; y_order <= max_order-x_order; y_order++) {
266 
267  // evaluate exact solution
268  FieldContainer<double> exact_solution(1, numInterpPoints);
269  u_exact(exact_solution, interp_points, x_order, y_order);
270 
271  int total_order = std::max(x_order + y_order, 1);
272 
273  for (int basis_order=total_order; basis_order <= max_order; basis_order++) {
274 
275  // set test tolerance
276  double zero = basis_order*basis_order*100*INTREPID_TOL;
277 
278  //create basis
279  Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
280  Teuchos::rcp(new Basis_HGRAD_TRI_Cn_FEM<double,FieldContainer<double> >(basis_order, pointtype[ptype]) );
281  int numFields = basis->getCardinality();
282 
283  // create cubatures
284  Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
285  Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*basis_order);
286  int numCubPointsCell = cellCub->getNumPoints();
287  int numCubPointsSide = sideCub->getNumPoints();
288 
289  /* Computational arrays. */
290  /* Section 1: Related to parent cell integration. */
291  FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
292  FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
293  FieldContainer<double> cub_weights_cell(numCubPointsCell);
294  FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
295  FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
296  FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
297  FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
298 
299  FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
300  FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
301  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
302  FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
303  FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
304  FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
305  FieldContainer<double> fe_matrix(1, numFields, numFields);
306 
307  FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
308  FieldContainer<double> rhs_and_soln_vector(1, numFields);
309 
310  /* Section 2: Related to subcell (side) integration. */
311  unsigned numSides = 3;
312  FieldContainer<double> cub_points_side(numCubPointsSide, sideDim);
313  FieldContainer<double> cub_weights_side(numCubPointsSide);
314  FieldContainer<double> cub_points_side_refcell(numCubPointsSide, cellDim);
315  FieldContainer<double> cub_points_side_physical(1, numCubPointsSide, cellDim);
316  FieldContainer<double> jacobian_side_refcell(1, numCubPointsSide, cellDim, cellDim);
317  FieldContainer<double> jacobian_det_side_refcell(1, numCubPointsSide);
318  FieldContainer<double> weighted_measure_side_refcell(1, numCubPointsSide);
319 
320  FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields, numCubPointsSide);
321  FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
322  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
323  FieldContainer<double> neumann_data_at_cub_points_side_physical(1, numCubPointsSide);
324  FieldContainer<double> neumann_fields_per_side(1, numFields);
325 
326  /* Section 3: Related to global interpolant. */
327  FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
328  FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
329  FieldContainer<double> interpolant(1, numInterpPoints);
330 
331  FieldContainer<int> ipiv(numFields);
332 
333 
334 
335  /******************* START COMPUTATION ***********************/
336 
337  // get cubature points and weights
338  cellCub->getCubature(cub_points_cell, cub_weights_cell);
339 
340  // compute geometric cell information
341  CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes, cell);
342  CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
343  CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
344 
345  // compute weighted measure
346  FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
347 
349  // Computing mass matrices:
350  // tabulate values of basis functions at (reference) cubature points
351  basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
352 
353  // transform values of basis functions
354  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
355  value_of_basis_at_cub_points_cell);
356 
357  // multiply with weighted measure
358  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
359  weighted_measure_cell,
360  transformed_value_of_basis_at_cub_points_cell);
361 
362  // compute mass matrices
363  FunctionSpaceTools::integrate<double>(fe_matrix,
364  transformed_value_of_basis_at_cub_points_cell,
365  weighted_transformed_value_of_basis_at_cub_points_cell,
366  COMP_BLAS);
368 
370  // Computing stiffness matrices:
371  // tabulate gradients of basis functions at (reference) cubature points
372  basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
373 
374  // transform gradients of basis functions
375  FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
376  jacobian_inv_cell,
377  grad_of_basis_at_cub_points_cell);
378 
379  // multiply with weighted measure
380  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
381  weighted_measure_cell,
382  transformed_grad_of_basis_at_cub_points_cell);
383 
384  // compute stiffness matrices and sum into fe_matrix
385  FunctionSpaceTools::integrate<double>(fe_matrix,
386  transformed_grad_of_basis_at_cub_points_cell,
387  weighted_transformed_grad_of_basis_at_cub_points_cell,
388  COMP_BLAS,
389  true);
391 
393  // Computing RHS contributions:
394  // map cell (reference) cubature points to physical space
395  CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes, cell);
396 
397  // evaluate rhs function
398  rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order);
399 
400  // compute rhs
401  FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
402  rhs_at_cub_points_cell_physical,
403  weighted_transformed_value_of_basis_at_cub_points_cell,
404  COMP_BLAS);
405 
406  // compute neumann b.c. contributions and adjust rhs
407  sideCub->getCubature(cub_points_side, cub_weights_side);
408  for (unsigned i=0; i<numSides; i++) {
409  // compute geometric cell information
410  CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
411  CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes, cell);
412  CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
413 
414  // compute weighted edge measure
415  FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
416  jacobian_side_refcell,
417  cub_weights_side,
418  i,
419  cell);
420 
421  // tabulate values of basis functions at side cubature points, in the reference parent cell domain
422  basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
423  // transform
424  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
425  value_of_basis_at_cub_points_side_refcell);
426 
427  // multiply with weighted measure
428  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
429  weighted_measure_side_refcell,
430  transformed_value_of_basis_at_cub_points_side_refcell);
431 
432  // compute Neumann data
433  // map side cubature points in reference parent cell domain to physical space
434  CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes, cell);
435  // now compute data
436  neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
437  cell, (int)i, x_order, y_order);
438 
439  FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
440  neumann_data_at_cub_points_side_physical,
441  weighted_transformed_value_of_basis_at_cub_points_side_refcell,
442  COMP_BLAS);
443 
444  // adjust RHS
445  RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
446  }
448 
450  // Solution of linear system:
451  int info = 0;
452  Teuchos::LAPACK<int, double> solver;
453  solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
455 
457  // Building interpolant:
458  // evaluate basis at interpolation points
459  basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
460  // transform values of basis functions
461  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
462  value_of_basis_at_interp_points);
463  FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
465 
466  /******************* END COMPUTATION ***********************/
467 
468  RealSpaceTools<double>::subtract(interpolant, exact_solution);
469 
470  *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
471  << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
472  << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
473  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
474 
475  if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
476  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
477  *outStream << "\n\nPatch test failed for solution polynomial order ("
478  << x_order << ", " << y_order << ") and basis order " << basis_order << "\n\n";
479  errorFlag++;
480  }
481  } // end for basis_order
482  } // end for y_order
483  } // end for x_order
484  } // end for ptype
485 
486  }
487  // Catch unexpected errors
488  catch (const std::logic_error & err) {
489  *outStream << err.what() << "\n\n";
490  errorFlag = -1000;
491  };
492 
493  if (errorFlag != 0)
494  std::cout << "End Result: TEST FAILED\n";
495  else
496  std::cout << "End Result: TEST PASSED\n";
497 
498  // reset format state of std::cout
499  std::cout.copyfmt(oldFormatState);
500 
501  return errorFlag;
502 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for the Intrepid::HGRAD_TRI_Cn_FEM class.
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.
Header file for the Intrepid::FunctionSpaceTools class.
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell...
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: