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