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