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