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