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, 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_HEX_I2_FEM) |\n" \
214  << "| |\n" \
215  << "| 1) Patch test involving mass and stiffness matrices, |\n" \
216  << "| for the Neumann problem on a physical parallelepiped |\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  << "| |\n" \
222  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
223  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
224  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
225  << "| |\n" \
226  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
227  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
228  << "| |\n" \
229  << "===============================================================================\n"\
230  << "| TEST 1: Patch test |\n"\
231  << "===============================================================================\n";
232 
233 
234  int errorFlag = 0;
235 
236  outStream -> precision(16);
237 
238 
239  try {
240 
241  int max_order = 2; // max total order of polynomial solution
242  DefaultCubatureFactory<double> cubFactory; // create factory
243  shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<> >()); // create parent cell topology
244  shards::CellTopology side(shards::getCellTopologyData< shards::Quadrilateral<> >()); // create relevant subcell (side) topology
245  int cellDim = cell.getDimension();
246  int sideDim = side.getDimension();
247  unsigned numSides = 6;
248 
249  // Define array containing points at which the solution is evaluated, on the reference tet.
250  int numIntervals = 10;
251  int numInterpPoints = (numIntervals + 1)*(numIntervals + 1)*(numIntervals + 1);
252  FieldContainer<double> interp_points_ref(numInterpPoints, 3);
253  int counter = 0;
254  for (int k=0; k<=numIntervals; k++) {
255  for (int j=0; j<=numIntervals; j++) {
256  for (int i=0; i<=numIntervals; i++) {
257  interp_points_ref(counter,0) = i*(1.0/numIntervals)-1.0;
258  interp_points_ref(counter,1) = j*(1.0/numIntervals)-1.0;
259  interp_points_ref(counter,2) = k*(1.0/numIntervals)-1.0;
260  counter++;
261  }
262  }
263  }
264 
265  /* Parent cell definition. */
266  FieldContainer<double> cell_nodes(1, 8, cellDim);
267 
268  // Generic parallelepiped.
269  cell_nodes(0, 0, 0) = -5.0;
270  cell_nodes(0, 0, 1) = -1.0;
271  cell_nodes(0, 0, 2) = 0.0;
272  cell_nodes(0, 1, 0) = 4.0;
273  cell_nodes(0, 1, 1) = 1.0;
274  cell_nodes(0, 1, 2) = 1.0;
275  cell_nodes(0, 2, 0) = 8.0;
276  cell_nodes(0, 2, 1) = 3.0;
277  cell_nodes(0, 2, 2) = 1.0;
278  cell_nodes(0, 3, 0) = -1.0;
279  cell_nodes(0, 3, 1) = 1.0;
280  cell_nodes(0, 3, 2) = 0.0;
281  cell_nodes(0, 4, 0) = 5.0;
282  cell_nodes(0, 4, 1) = 9.0;
283  cell_nodes(0, 4, 2) = 1.0;
284  cell_nodes(0, 5, 0) = 14.0;
285  cell_nodes(0, 5, 1) = 11.0;
286  cell_nodes(0, 5, 2) = 2.0;
287  cell_nodes(0, 6, 0) = 18.0;
288  cell_nodes(0, 6, 1) = 13.0;
289  cell_nodes(0, 6, 2) = 2.0;
290  cell_nodes(0, 7, 0) = 9.0;
291  cell_nodes(0, 7, 1) = 11.0;
292  cell_nodes(0, 7, 2) = 1.0;
293 
294  // Reference hex.
295  /*cell_nodes(0, 0, 0) = -1.0;
296  cell_nodes(0, 0, 1) = -1.0;
297  cell_nodes(0, 0, 2) = -1.0;
298  cell_nodes(0, 1, 0) = 1.0;
299  cell_nodes(0, 1, 1) = -1.0;
300  cell_nodes(0, 1, 2) = -1.0;
301  cell_nodes(0, 2, 0) = 1.0;
302  cell_nodes(0, 2, 1) = 1.0;
303  cell_nodes(0, 2, 2) = -1.0;
304  cell_nodes(0, 3, 0) = -1.0;
305  cell_nodes(0, 3, 1) = 1.0;
306  cell_nodes(0, 3, 2) = -1.0;
307  cell_nodes(0, 4, 0) = -1.0;
308  cell_nodes(0, 4, 1) = -1.0;
309  cell_nodes(0, 4, 2) = 1.0;
310  cell_nodes(0, 5, 0) = 1.0;
311  cell_nodes(0, 5, 1) = -1.0;
312  cell_nodes(0, 5, 2) = 1.0;
313  cell_nodes(0, 6, 0) = 1.0;
314  cell_nodes(0, 6, 1) = 1.0;
315  cell_nodes(0, 6, 2) = 1.0;
316  cell_nodes(0, 7, 0) = -1.0;
317  cell_nodes(0, 7, 1) = 1.0;
318  cell_nodes(0, 7, 2) = 1.0;*/
319 
320  FieldContainer<double> interp_points(1, numInterpPoints, cellDim);
321  CellTools<double>::mapToPhysicalFrame(interp_points, interp_points_ref, cell_nodes, cell);
322  interp_points.resize(numInterpPoints, cellDim);
323 
324  for (int x_order=0; x_order <= max_order; x_order++) {
325  for (int y_order=0; y_order <= max_order - x_order; y_order++) {
326  for (int z_order=0; z_order <= max_order - x_order - y_order; z_order++) {
327 
328  // evaluate exact solution
329  FieldContainer<double> exact_solution(1, numInterpPoints);
330  u_exact(exact_solution, interp_points, x_order, y_order, z_order);
331 
332  int basis_order = 2;
333 
334  // set test tolerance;
335  double zero = basis_order*basis_order*basis_order*100*INTREPID_TOL;
336 
337  //create basis
338  Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
339  Teuchos::rcp(new Basis_HGRAD_HEX_I2_FEM<double,FieldContainer<double> >() );
340  int numFields = basis->getCardinality();
341 
342  // create cubatures
343  Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
344  Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*basis_order);
345  int numCubPointsCell = cellCub->getNumPoints();
346  int numCubPointsSide = sideCub->getNumPoints();
347 
348  /* Computational arrays. */
349  /* Section 1: Related to parent cell integration. */
350  FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
351  FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
352  FieldContainer<double> cub_weights_cell(numCubPointsCell);
353  FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
354  FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
355  FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
356  FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
357 
358  FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
359  FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
360  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
361  FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
362  FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
363  FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
364  FieldContainer<double> fe_matrix(1, numFields, numFields);
365 
366  FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
367  FieldContainer<double> rhs_and_soln_vector(1, numFields);
368 
369  /* Section 2: Related to subcell (side) integration. */
370  FieldContainer<double> cub_points_side(numCubPointsSide, sideDim);
371  FieldContainer<double> cub_weights_side(numCubPointsSide);
372  FieldContainer<double> cub_points_side_refcell(numCubPointsSide, cellDim);
373  FieldContainer<double> cub_points_side_physical(1, numCubPointsSide, cellDim);
374  FieldContainer<double> jacobian_side_refcell(1, numCubPointsSide, cellDim, cellDim);
375  FieldContainer<double> jacobian_det_side_refcell(1, numCubPointsSide);
376  FieldContainer<double> weighted_measure_side_refcell(1, numCubPointsSide);
377 
378  FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields, numCubPointsSide);
379  FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
380  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
381  FieldContainer<double> neumann_data_at_cub_points_side_physical(1, numCubPointsSide);
382  FieldContainer<double> neumann_fields_per_side(1, numFields);
383 
384  /* Section 3: Related to global interpolant. */
385  FieldContainer<double> value_of_basis_at_interp_points_ref(numFields, numInterpPoints);
386  FieldContainer<double> transformed_value_of_basis_at_interp_points_ref(1, numFields, numInterpPoints);
387  FieldContainer<double> interpolant(1, numInterpPoints);
388 
389  FieldContainer<int> ipiv(numFields);
390 
391 
392 
393  /******************* START COMPUTATION ***********************/
394 
395  // get cubature points and weights
396  cellCub->getCubature(cub_points_cell, cub_weights_cell);
397 
398  // compute geometric cell information
399  CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes, cell);
400  CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
401  CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
402 
403  // compute weighted measure
404  FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
405 
407  // Computing mass matrices:
408  // tabulate values of basis functions at (reference) cubature points
409  basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
410 
411  // transform values of basis functions
412  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
413  value_of_basis_at_cub_points_cell);
414 
415  // multiply with weighted measure
416  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
417  weighted_measure_cell,
418  transformed_value_of_basis_at_cub_points_cell);
419 
420  // compute mass matrices
421  FunctionSpaceTools::integrate<double>(fe_matrix,
422  transformed_value_of_basis_at_cub_points_cell,
423  weighted_transformed_value_of_basis_at_cub_points_cell,
424  COMP_BLAS);
426 
428  // Computing stiffness matrices:
429  // tabulate gradients of basis functions at (reference) cubature points
430  basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
431 
432  // transform gradients of basis functions
433  FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
434  jacobian_inv_cell,
435  grad_of_basis_at_cub_points_cell);
436 
437  // multiply with weighted measure
438  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
439  weighted_measure_cell,
440  transformed_grad_of_basis_at_cub_points_cell);
441 
442  // compute stiffness matrices and sum into fe_matrix
443  FunctionSpaceTools::integrate<double>(fe_matrix,
444  transformed_grad_of_basis_at_cub_points_cell,
445  weighted_transformed_grad_of_basis_at_cub_points_cell,
446  COMP_BLAS,
447  true);
449 
451  // Computing RHS contributions:
452  // map cell (reference) cubature points to physical space
453  CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes, cell);
454 
455  // evaluate rhs function
456  rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order, z_order);
457 
458  // compute rhs
459  FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
460  rhs_at_cub_points_cell_physical,
461  weighted_transformed_value_of_basis_at_cub_points_cell,
462  COMP_BLAS);
463 
464  // compute neumann b.c. contributions and adjust rhs
465  sideCub->getCubature(cub_points_side, cub_weights_side);
466  for (unsigned i=0; i<numSides; i++) {
467  // compute geometric cell information
468  CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
469  CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes, cell);
470  CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
471 
472  // compute weighted face measure
473  FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_side_refcell,
474  jacobian_side_refcell,
475  cub_weights_side,
476  i,
477  cell);
478 
479  // tabulate values of basis functions at side cubature points, in the reference parent cell domain
480  basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
481  // transform
482  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
483  value_of_basis_at_cub_points_side_refcell);
484 
485  // multiply with weighted measure
486  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
487  weighted_measure_side_refcell,
488  transformed_value_of_basis_at_cub_points_side_refcell);
489 
490  // compute Neumann data
491  // map side cubature points in reference parent cell domain to physical space
492  CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes, cell);
493  // now compute data
494  neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
495  cell, (int)i, x_order, y_order, z_order);
496 
497  FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
498  neumann_data_at_cub_points_side_physical,
499  weighted_transformed_value_of_basis_at_cub_points_side_refcell,
500  COMP_BLAS);
501 
502  // adjust RHS
503  RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
504  }
506 
508  // Solution of linear system:
509  int info = 0;
510  Teuchos::LAPACK<int, double> solver;
511  solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
513 
515  // Building interpolant:
516  // evaluate basis at interpolation points
517  basis->getValues(value_of_basis_at_interp_points_ref, interp_points_ref, OPERATOR_VALUE);
518  // transform values of basis functions
519  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points_ref,
520  value_of_basis_at_interp_points_ref);
521  FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points_ref);
523 
524  /******************* END COMPUTATION ***********************/
525 
526  RealSpaceTools<double>::subtract(interpolant, exact_solution);
527 
528  *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
529  << x_order << ", " << y_order << ", " << z_order
530  << ") and finite element interpolant of order " << basis_order << ": "
531  << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
532  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
533 
534  if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
535  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
536  *outStream << "\n\nPatch test failed for solution polynomial order ("
537  << x_order << ", " << y_order << ", " << z_order << ") and basis order " << basis_order << "\n\n";
538  errorFlag++;
539  }
540  } // end for z_order
541  } // end for y_order
542  } // end for x_order
543 
544  }
545  // Catch unexpected errors
546  catch (const std::logic_error & err) {
547  *outStream << err.what() << "\n\n";
548  errorFlag = -1000;
549  };
550 
551  if (errorFlag != 0)
552  std::cout << "End Result: TEST FAILED\n";
553  else
554  std::cout << "End Result: TEST PASSED\n";
555 
556  // reset format state of std::cout
557  std::cout.copyfmt(oldFormatState);
558 
559  return errorFlag;
560 }
Implementation of basic linear algebra functionality in Euclidean space.
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
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::HGRAD_HEX_I2_FEM class.
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: