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_PYR_I2_FEM) |\n" \
214  << "| |\n" \
215  << "| 1) Patch test involving mass and stiffness matrices, |\n" \
216  << "| for the Neumann problem on a pyramid 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  << "| Mauro Perego (mperego@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::Pyramid<> >()); // create parent cell topology
244  shards::CellTopology sideQ(shards::getCellTopologyData< shards::Quadrilateral<> >()); // create relevant subcell (side) topology
245  shards::CellTopology sideT(shards::getCellTopologyData< shards::Triangle<> >());
246  int cellDim = cell.getDimension();
247  int sideQDim = sideQ.getDimension();
248  int sideTDim = sideT.getDimension();
249 
250  // Define array containing points at which the solution is evaluated, on the reference Pyramid.
251  int numIntervals = 10;
252  int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals + 3))/6;
253  FieldContainer<double> interp_points_ref(numInterpPoints, 3);
254  int counter = 0;
255  for (int k=0; k<=numIntervals; k++) {
256  for (int j=0; j<=numIntervals; j++) {
257  for (int i=0; i<=numIntervals; i++) {
258  if (i+j+k <= numIntervals) {
259  interp_points_ref(counter,0) = i*(1.0/numIntervals);
260  interp_points_ref(counter,1) = j*(1.0/numIntervals);
261  interp_points_ref(counter,2) = k*(1.0/numIntervals);
262  counter++;
263  }
264  }
265  }
266  }
267 
268  /* Definition of parent cell. */
269  FieldContainer<double> cell_nodes(1, 5, cellDim);
270 
271  // Pyramid with affine mapping
272 
273  cell_nodes(0, 0, 0) = -4.0;
274  cell_nodes(0, 0, 1) = -9.0;
275  cell_nodes(0, 0, 2) = -5.0;
276  cell_nodes(0, 1, 0) = -6.0;
277  cell_nodes(0, 1, 1) = -3.0;
278  cell_nodes(0, 1, 2) = 3.0;
279  cell_nodes(0, 2, 0) = 10.0;
280  cell_nodes(0, 2, 1) = 5.0;
281  cell_nodes(0, 2, 2) = 7.0;
282  cell_nodes(0, 3, 0) = 12.0;
283  cell_nodes(0, 3, 1) = -1.0;
284  cell_nodes(0, 3, 2) = -1.0;
285  cell_nodes(0, 4, 0) = 5.0;
286  cell_nodes(0, 4, 1) = 5.0;
287  cell_nodes(0, 4, 2) = -3.0;
288 
289 
290 /*
291  cell_nodes(0, 0, 0) = 0.0;
292  cell_nodes(0, 0, 1) = -6.0;
293  cell_nodes(0, 0, 2) = -2.0;
294  cell_nodes(0, 1, 0) = 8.0;
295  cell_nodes(0, 1, 1) =-10.0;
296  cell_nodes(0, 1, 2) = 0.0;
297  cell_nodes(0, 2, 0) = 6.0;
298  cell_nodes(0, 2, 1) = 2.0;
299  cell_nodes(0, 2, 2) = 4.0;
300  cell_nodes(0, 3, 0) = -2.0;
301  cell_nodes(0, 3, 1) = 6.0;
302  cell_nodes(0, 3, 2) = 2.0;
303  cell_nodes(0, 4, 0) = 6.0;
304  cell_nodes(0, 4, 1) = 5.0;
305  cell_nodes(0, 4, 2) = -3.0;
306 */
307 
308  // reference Pyramid
309  /*cell_nodes(0, 0, 0) = -1.0;
310  cell_nodes(0, 0, 1) = -1.0;
311  cell_nodes(0, 0, 2) = 0.0;
312  cell_nodes(0, 1, 0) = 1.0;
313  cell_nodes(0, 1, 1) = -1.0;
314  cell_nodes(0, 1, 2) = 0.0;
315  cell_nodes(0, 2, 0) = 1.0;
316  cell_nodes(0, 2, 1) = 1.0;
317  cell_nodes(0, 2, 2) = 0.0;
318  cell_nodes(0, 3, 0) = -1.0;
319  cell_nodes(0, 3, 1) = 1.0;
320  cell_nodes(0, 3, 2) = 0.0;
321  cell_nodes(0, 4, 0) = 0.0;
322  cell_nodes(0, 4, 1) = 0.0;
323  cell_nodes(0, 4, 2) = 1.0;*/
324 
325 
326  FieldContainer<double> interp_points(1, numInterpPoints, cellDim);
327  CellTools<double>::mapToPhysicalFrame(interp_points, interp_points_ref, cell_nodes, cell);
328  interp_points.resize(numInterpPoints, cellDim);
329 
330  for (int x_order=0; x_order <= max_order; x_order++) {
331  for (int y_order=0; y_order <= max_order-x_order; y_order++) {
332  for (int z_order=0; z_order <= max_order-x_order-y_order; z_order++) {
333 
334  // evaluate exact solution
335  FieldContainer<double> exact_solution(1, numInterpPoints);
336  u_exact(exact_solution, interp_points, x_order, y_order, z_order);
337 
338  int basis_order = 2;
339 
340  // set test tolerance;
341  double zero = basis_order*basis_order*basis_order*100*INTREPID_TOL;
342 
343  //create basis
344  Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
345  Teuchos::rcp(new Basis_HGRAD_PYR_I2_FEM<double,FieldContainer<double> >() );
346  int numFields = basis->getCardinality();
347 
348  // create cubatures
349  Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
350  Teuchos::RCP<Cubature<double> > sideQCub = cubFactory.create(sideQ, 2*basis_order);
351  Teuchos::RCP<Cubature<double> > sideTCub = cubFactory.create(sideT, 2*basis_order);
352  int numCubPointsCell = cellCub->getNumPoints();
353  int numCubPointsSideQ = sideQCub->getNumPoints();
354  int numCubPointsSideT = sideTCub->getNumPoints();
355 
356  /* Computational arrays. */
357  /* Section 1: Related to parent cell integration. */
358  FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
359  FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
360  FieldContainer<double> cub_weights_cell(numCubPointsCell);
361  FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
362  FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
363  FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
364  FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
365 
366  FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
367  FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
368  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
369  FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
370  FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
371  FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
372  FieldContainer<double> fe_matrix(1, numFields, numFields);
373 
374  FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
375  FieldContainer<double> rhs_and_soln_vector(1, numFields);
376 
377  /* Section 2: Related to subcell (side) integration. */
378  unsigned numSides = 5;
379  unsigned numSidesT = 4;
380  FieldContainer<double> cub_points_sideQ(numCubPointsSideQ, sideQDim);
381  FieldContainer<double> cub_points_sideT(numCubPointsSideT, sideTDim);
382  FieldContainer<double> cub_weights_sideQ(numCubPointsSideQ);
383  FieldContainer<double> cub_weights_sideT(numCubPointsSideT);
384  FieldContainer<double> cub_points_sideQ_refcell(numCubPointsSideQ, cellDim);
385  FieldContainer<double> cub_points_sideT_refcell(numCubPointsSideT, cellDim);
386  FieldContainer<double> cub_points_sideQ_physical(1, numCubPointsSideQ, cellDim);
387  FieldContainer<double> cub_points_sideT_physical(1, numCubPointsSideT, cellDim);
388  FieldContainer<double> jacobian_sideQ_refcell(1, numCubPointsSideQ, cellDim, cellDim);
389  FieldContainer<double> jacobian_sideT_refcell(1, numCubPointsSideT, cellDim, cellDim);
390  FieldContainer<double> jacobian_det_sideQ_refcell(1, numCubPointsSideQ);
391  FieldContainer<double> jacobian_det_sideT_refcell(1, numCubPointsSideT);
392  FieldContainer<double> weighted_measure_sideQ_refcell(1, numCubPointsSideQ);
393  FieldContainer<double> weighted_measure_sideT_refcell(1, numCubPointsSideT);
394 
395  FieldContainer<double> value_of_basis_at_cub_points_sideQ_refcell(numFields, numCubPointsSideQ);
396  FieldContainer<double> value_of_basis_at_cub_points_sideT_refcell(numFields, numCubPointsSideT);
397  FieldContainer<double> transformed_value_of_basis_at_cub_points_sideQ_refcell(1, numFields, numCubPointsSideQ);
398  FieldContainer<double> transformed_value_of_basis_at_cub_points_sideT_refcell(1, numFields, numCubPointsSideT);
399  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_sideQ_refcell(1, numFields, numCubPointsSideQ);
400  FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_sideT_refcell(1, numFields, numCubPointsSideT);
401  FieldContainer<double> neumann_data_at_cub_points_sideQ_physical(1, numCubPointsSideQ);
402  FieldContainer<double> neumann_data_at_cub_points_sideT_physical(1, numCubPointsSideT);
403  FieldContainer<double> neumann_fields_per_side(1, numFields);
404 
405  /* Section 3: Related to global interpolant. */
406  FieldContainer<double> value_of_basis_at_interp_points_ref(numFields, numInterpPoints);
407  FieldContainer<double> transformed_value_of_basis_at_interp_points_ref(1, numFields, numInterpPoints);
408  FieldContainer<double> interpolant(1, numInterpPoints);
409 
410  FieldContainer<int> ipiv(numFields);
411 
412 
413 
414  /******************* START COMPUTATION ***********************/
415 
416  // get cubature points and weights
417  cellCub->getCubature(cub_points_cell, cub_weights_cell);
418 
419  // compute geometric cell information
420  CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes, cell);
421  CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
422  CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
423 
424  // compute weighted measure
425  FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
426 
428  // Computing mass matrices:
429  // tabulate values of basis functions at (reference) cubature points
430  basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
431 
432  // transform values of basis functions
433  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
434  value_of_basis_at_cub_points_cell);
435 
436  // multiply with weighted measure
437  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
438  weighted_measure_cell,
439  transformed_value_of_basis_at_cub_points_cell);
440 
441  // compute mass matrices
442  FunctionSpaceTools::integrate<double>(fe_matrix,
443  transformed_value_of_basis_at_cub_points_cell,
444  weighted_transformed_value_of_basis_at_cub_points_cell,
445  COMP_BLAS);
447 
449  // Computing stiffness matrices:
450  // tabulate gradients of basis functions at (reference) cubature points
451  basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
452 
453  // transform gradients of basis functions
454  FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
455  jacobian_inv_cell,
456  grad_of_basis_at_cub_points_cell);
457 
458  // multiply with weighted measure
459  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
460  weighted_measure_cell,
461  transformed_grad_of_basis_at_cub_points_cell);
462 
463  // compute stiffness matrices and sum into fe_matrix
464  FunctionSpaceTools::integrate<double>(fe_matrix,
465  transformed_grad_of_basis_at_cub_points_cell,
466  weighted_transformed_grad_of_basis_at_cub_points_cell,
467  COMP_BLAS,
468  true);
470 
472  // Computing RHS contributions:
473  // map cell (reference) cubature points to physical space
474  CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes, cell);
475 
476  // evaluate rhs function
477  rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order, z_order);
478 
479  // compute rhs
480  FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
481  rhs_at_cub_points_cell_physical,
482  weighted_transformed_value_of_basis_at_cub_points_cell,
483  COMP_BLAS);
484 
485  // compute neumann b.c. contributions and adjust rhs
486  sideQCub->getCubature(cub_points_sideQ, cub_weights_sideQ);
487  sideTCub->getCubature(cub_points_sideT, cub_weights_sideT);
488 
489  for (unsigned i=0; i<numSidesT; i++) {
490  // compute geometric cell information
491  CellTools<double>::mapToReferenceSubcell(cub_points_sideT_refcell, cub_points_sideT, sideTDim, (int)i, cell);
492  CellTools<double>::setJacobian(jacobian_sideT_refcell, cub_points_sideT_refcell, cell_nodes, cell);
493  CellTools<double>::setJacobianDet(jacobian_det_sideT_refcell, jacobian_sideT_refcell);
494 
495  // compute weighted face measure
496  FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_sideT_refcell,
497  jacobian_sideT_refcell,
498  cub_weights_sideT,
499  i,
500  cell);
501 
502  // tabulate values of basis functions at side cubature points, in the reference parent cell domain
503  basis->getValues(value_of_basis_at_cub_points_sideT_refcell, cub_points_sideT_refcell, OPERATOR_VALUE);
504  // transform
505  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_sideT_refcell,
506  value_of_basis_at_cub_points_sideT_refcell);
507 
508  // multiply with weighted measure
509  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_sideT_refcell,
510  weighted_measure_sideT_refcell,
511  transformed_value_of_basis_at_cub_points_sideT_refcell);
512 
513  // compute Neumann data
514  // map side cubature points in reference parent cell domain to physical space
515  CellTools<double>::mapToPhysicalFrame(cub_points_sideT_physical, cub_points_sideT_refcell, cell_nodes, cell);
516  // now compute data
517  neumann(neumann_data_at_cub_points_sideT_physical, cub_points_sideT_physical, jacobian_sideT_refcell,
518  cell, (int)i, x_order, y_order, z_order);
519 
520  FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
521  neumann_data_at_cub_points_sideT_physical,
522  weighted_transformed_value_of_basis_at_cub_points_sideT_refcell,
523  COMP_BLAS);
524 
525  // adjust RHS
526  RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);
527  }
528 
529  for (unsigned i=numSidesT; i<numSides; i++) {
530  // compute geometric cell information
531  CellTools<double>::mapToReferenceSubcell(cub_points_sideQ_refcell, cub_points_sideQ, sideQDim, (int)i, cell);
532  CellTools<double>::setJacobian(jacobian_sideQ_refcell, cub_points_sideQ_refcell, cell_nodes, cell);
533  CellTools<double>::setJacobianDet(jacobian_det_sideQ_refcell, jacobian_sideQ_refcell);
534 
535  // compute weighted face measure
536  FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_sideQ_refcell,
537  jacobian_sideQ_refcell,
538  cub_weights_sideQ,
539  i,
540  cell);
541 
542  // tabulate values of basis functions at side cubature points, in the reference parent cell domain
543  basis->getValues(value_of_basis_at_cub_points_sideQ_refcell, cub_points_sideQ_refcell, OPERATOR_VALUE);
544  // transform
545  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_sideQ_refcell,
546  value_of_basis_at_cub_points_sideQ_refcell);
547 
548  // multiply with weighted measure
549  FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_sideQ_refcell,
550  weighted_measure_sideQ_refcell,
551  transformed_value_of_basis_at_cub_points_sideQ_refcell);
552 
553  // compute Neumann data
554  // map side cubature points in reference parent cell domain to physical space
555  CellTools<double>::mapToPhysicalFrame(cub_points_sideQ_physical, cub_points_sideQ_refcell, cell_nodes, cell);
556  // now compute data
557  neumann(neumann_data_at_cub_points_sideQ_physical, cub_points_sideQ_physical, jacobian_sideQ_refcell,
558  cell, (int)i, x_order, y_order, z_order);
559 
560  FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
561  neumann_data_at_cub_points_sideQ_physical,
562  weighted_transformed_value_of_basis_at_cub_points_sideQ_refcell,
563  COMP_BLAS);
564 
565  // adjust RHS
566  RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
567  }
569 
571  // Solution of linear system:
572  int info = 0;
573  Teuchos::LAPACK<int, double> solver;
574  solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
576 
577 // std::cout << rhs_and_soln_vector;
578 
580  // Building interpolant:
581  // evaluate basis at interpolation points
582  basis->getValues(value_of_basis_at_interp_points_ref, interp_points_ref, OPERATOR_VALUE);
583  // transform values of basis functions
584  FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points_ref,
585  value_of_basis_at_interp_points_ref);
586  FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points_ref);
588 
589  /******************* END COMPUTATION ***********************/
590 
591  RealSpaceTools<double>::subtract(interpolant, exact_solution);
592 
593  *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
594  << x_order << ", " << y_order << ", " << z_order
595  << ") and finite element interpolant of order " << basis_order << ": "
596  << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
597  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
598 
599  if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
600  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
601  *outStream << "\n\nPatch test failed for solution polynomial order ("
602  << x_order << ", " << y_order << ", " << z_order << ") and basis order " << basis_order << "\n\n";
603  errorFlag++;
604  }
605  } // end for z_order
606  } // end for y_order
607  } // end for x_order
608 
609  }
610  // Catch unexpected errors
611  catch (const std::logic_error & err) {
612  *outStream << err.what() << "\n\n";
613  errorFlag = -1000;
614  };
615 
616  if (errorFlag != 0)
617  std::cout << "End Result: TEST FAILED\n";
618  else
619  std::cout << "End Result: TEST PASSED\n";
620 
621  // reset format state of std::cout
622  std::cout.copyfmt(oldFormatState);
623 
624  return errorFlag;
625 }
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...
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
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.
Header file for the Intrepid::HGRAD_PYR_I2_FEM class.
A stateless class for operations on cell data. Provides methods for: