Intrepid
test_03.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
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
63using namespace std;
64using namespace Intrepid;
65
66void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int, int );
67void bcFunc( FieldContainer<double> &, const FieldContainer<double> & ,
68 const shards::CellTopology & ,
69 int , int , int , int , int );
70void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int );
71
72void u_exact( FieldContainer<double> &result,
73 const FieldContainer<double> &points,
74 int comp,
75 int xd,
76 int yd,
77 int zd)
78{
79 for (int cell=0;cell<result.dimension(0);cell++){
80 for (int pt=0;pt<result.dimension(1);pt++) {
81 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
82 *std::pow(points(cell,pt,2),zd);
83 }
84 }
85 return;
86}
87
88void bcFunc( FieldContainer<double> & result,
89 const FieldContainer<double> & points ,
90 const shards::CellTopology & parentCell ,
91 int sideOrdinal , int comp , int a, int b, int c )
92{
93
94 int numCells = result.dimension(0);
95 int numPoints = result.dimension(1);
96
97 // reference face normal
98 FieldContainer<double> normal(3);
99 CellTools<double>::getReferenceFaceNormal(normal,sideOrdinal,parentCell);
100
101 result.initialize();
102
103 if (comp == 0) {
104 for (int cell=0;cell<numCells;cell++) {
105 for (int pt=0;pt<numPoints;pt++) {
106 // first component
107 if (c > 0) {
108 result(cell,pt,0) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
109 }
110 if (b > 0) {
111 result(cell,pt,0) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
112 }
113 // second component
114 if (b > 0) {
115 result(cell,pt,1) = b * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
116 }
117 // third component
118 if (c > 0) {
119 result(cell,pt,2) = c * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
120 }
121 }
122 }
123 }
124 else if (comp == 1) {
125 for (int cell=0;cell<numCells;cell++) {
126 for (int pt=0;pt<numPoints;pt++) {
127 // first component
128 if (a > 0) {
129 result(cell,pt,0) = a * normal(1) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
130 }
131 // second component
132 if (c > 0) {
133 result(cell,pt,1) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
134 }
135 if (a > 0) {
136 result(cell,pt,1) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
137 }
138 // third component
139 if (c > 0) {
140 result(cell,pt,2) = c * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
141 }
142 }
143 }
144 }
145 else if (comp == 2) {
146 for (int cell=0;cell<numCells;cell++) {
147 for (int pt=0;pt<numPoints;pt++) {
148 // first component
149 if (a > 0) {
150 result(cell,pt,0) = a * normal(2) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
151 }
152 // second component
153 if (b > 0) {
154 result(cell,pt,1) = b * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
155 }
156 // third component
157 if (b > 0) {
158 result(cell,pt,2) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1)* std::pow(points(cell,pt,2),c);
159 }
160 if (a > 0) {
161 result(cell,pt,2) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
162 }
163 }
164 }
165 }
166
167 return;
168}
169
170void rhsFunc( FieldContainer<double> & result ,
171 const FieldContainer<double> &points ,
172 int comp,
173 int a,
174 int b,
175 int c )
176{
177 // rhs is curl(curl(E))+E, so load up the exact solution
178
179 u_exact(result,points,comp,a,b,c);
180
181 if (comp == 0) {
182 // component 0
183 if (b>=2) {
184 for (int cell=0;cell<result.dimension(0);cell++) {
185 for (int pt=0;pt<result.dimension(1);pt++) {
186 result(cell,pt,0) -= b*(b-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b-2)*std::pow(points(cell,pt,2),c);
187 }
188 }
189 }
190 if (c>=2) {
191 for (int cell=0;cell<result.dimension(0);cell++) {
192 for (int pt=0;pt<result.dimension(1);pt++) {
193 result(cell,pt,0)-=c*(c-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
194 }
195 }
196 }
197 // component one
198 if ( (a>=1) && (b>=1) ) {
199 for (int cell=0;cell<result.dimension(0);cell++) {
200 for (int pt=0;pt<result.dimension(1);pt++) {
201 result(cell,pt,1)+=a*b*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b-1)*std::pow(points(cell,pt,2),c);
202 }
203 }
204 }
205 // component two
206 if ( (a>=1) && (c>=1) ) {
207 for (int cell=0;cell<result.dimension(0);cell++) {
208 for (int pt=0;pt<result.dimension(1);pt++) {
209 result(cell,pt,2)+=a*c*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-1);
210 }
211 }
212 }
213 }
214 if (comp == 1) {
215 for (int cell=0;cell<result.dimension(0);cell++) {
216 for (int pt=0;pt<result.dimension(1);pt++) {
217 // first component
218 if (a > 0 && b > 0) {
219 result(cell,pt,0) += a * b * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
220 }
221 // second component
222 if (c > 1) {
223 result(cell,pt,1) -= c*(c-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
224 }
225 if (a > 1) {
226 result(cell,pt,1) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
227 }
228 // third component
229 if (b > 0 && c > 0) {
230 result(cell,pt,2) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
231 }
232 }
233 }
234 }
235 else if (comp == 2) {
236 for (int cell=0;cell<result.dimension(0);cell++) {
237 for (int pt=0;pt<result.dimension(1);pt++) {
238 // first component
239 if (a>0 && c>0) {
240 result(cell,pt,0) += a * c * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
241 }
242 // second component
243 if (b>0 && c>0) {
244 result(cell,pt,1) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
245 }
246 // third component
247 if (b>1) {
248 result(cell,pt,2) -= b*(b-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-2) * std::pow(points(cell,pt,2),c);
249 }
250 if (a>1) {
251 result(cell,pt,2) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
252 }
253 }
254 }
255 }
256 return;
257}
258
259int main(int argc, char *argv[]) {
260 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
261
262 // This little trick lets us print to std::cout only if
263 // a (dummy) command-line argument is provided.
264 int iprint = argc - 1;
265 Teuchos::RCP<std::ostream> outStream;
266 Teuchos::oblackholestream bhs; // outputs nothing
267 if (iprint > 0)
268 outStream = Teuchos::rcp(&std::cout, false);
269 else
270 outStream = Teuchos::rcp(&bhs, false);
271
272 // Save the format state of the original std::cout.
273 Teuchos::oblackholestream oldFormatState;
274 oldFormatState.copyfmt(std::cout);
275
276 *outStream \
277 << "===============================================================================\n" \
278 << "| |\n" \
279 << "| Unit Test (Basis_HGRAD_TRI_In_FEM) |\n" \
280 << "| |\n" \
281 << "| 1) Patch test involving H(curl) matrices |\n" \
282 << "| |\n" \
283 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
284 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
285 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
286 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
287 << "| |\n" \
288 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
289 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
290 << "| |\n" \
291 << "===============================================================================\n" \
292 << "| TEST 1: Patch test for mass + curl-curl matrices |\n" \
293 << "===============================================================================\n";
294
295
296 int errorFlag = 0;
297
298 outStream -> precision(16);
299
300 try {
301 DefaultCubatureFactory<double> cubFactory; // create cubature factory
302 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >()); // create parent cell topology
303
304 shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());
305
306 int cellDim = cell.getDimension();
307 int sideDim = side.getDimension();
308
309 int min_order = 1;
310 int max_order = 4;
311
312 int numIntervals = max_order;
313 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
314 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
315 int counter = 0;
316 for (int j=0; j<=numIntervals; j++) {
317 for (int i=0; i<=numIntervals-j; i++) {
318 for (int k=0;k<numIntervals-j-i;k++) {
319 interp_points_ref(counter,0) = i*(1.0/numIntervals);
320 interp_points_ref(counter,1) = j*(1.0/numIntervals);
321 interp_points_ref(counter,2) = k*(1.0/numIntervals);
322 counter++;
323 }
324 }
325 }
326
327 for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
328 // create basis
329 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
330 Teuchos::rcp(new Basis_HCURL_TET_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_EQUISPACED) );
331
332 int numFields = basis->getCardinality();
333
334 // create cubatures
335 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
336 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
337
338 int numCubPointsCell = cellCub->getNumPoints();
339 int numCubPointsSide = sideCub->getNumPoints();
340
341 // hold cubature information
342 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
343 FieldContainer<double> cub_weights_cell(numCubPointsCell);
344 FieldContainer<double> cub_points_side(numCubPointsSide,sideDim);
345 FieldContainer<double> cub_weights_side(numCubPointsSide);
346 FieldContainer<double> cub_points_side_refcell(numCubPointsSide,cellDim);
347
348 // hold basis function information on refcell
349 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
350 FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
351 FieldContainer<double> curl_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
352 FieldContainer<double> w_curl_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
353 FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields,numCubPointsSide,cellDim );
354 FieldContainer<double> w_value_of_basis_at_cub_points_side_refcell(1,numFields,numCubPointsSide,cellDim );
355
356
357 // holds rhs data
358 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim);
359 FieldContainer<double> bc_func_at_cub_points_side_refcell(1,numCubPointsSide,cellDim);
360 FieldContainer<double> bc_fields_per_side(1,numFields);
361
362 // FEM mass matrix
363 FieldContainer<double> fe_matrix_bak(1,numFields,numFields);
364 FieldContainer<double> fe_matrix(1,numFields,numFields);
365 FieldContainer<double> rhs_and_soln_vec(1,numFields);
366
367 FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim);
368 FieldContainer<double> interpolant( 1, numInterpPoints , cellDim );
369
370 int info = 0;
371 Teuchos::LAPACK<int, double> solver;
372
373 // set test tolerance
374 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
375
376 // build matrices outside the loop, and then just do the rhs
377 // for each iteration
378 cellCub->getCubature(cub_points_cell, cub_weights_cell);
379 sideCub->getCubature(cub_points_side, cub_weights_side);
380
381 // need the vector basis
382 basis->getValues(value_of_basis_at_cub_points_cell,
383 cub_points_cell,
384 OPERATOR_VALUE);
385 basis->getValues(curl_of_basis_at_cub_points_cell,
386 cub_points_cell,
387 OPERATOR_CURL);
388
389 basis->getValues( value_of_basis_at_interp_points ,
390 interp_points_ref ,
391 OPERATOR_VALUE );
392
393
394 // construct mass and curl-curl matrices
395 cub_weights_cell.resize(1,numCubPointsCell);
396 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
397 cub_weights_cell ,
398 value_of_basis_at_cub_points_cell );
399 FunctionSpaceTools::multiplyMeasure<double>(w_curl_of_basis_at_cub_points_cell ,
400 cub_weights_cell ,
401 curl_of_basis_at_cub_points_cell );
402 cub_weights_cell.resize(numCubPointsCell);
403
404
405 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
406 curl_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
407 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
408 w_value_of_basis_at_cub_points_cell ,
409 value_of_basis_at_cub_points_cell ,
410 COMP_BLAS );
411 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
412 w_curl_of_basis_at_cub_points_cell ,
413 curl_of_basis_at_cub_points_cell ,
414 COMP_BLAS ,
415 true );
416 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
417 curl_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
418
419 for (int comp=0;comp<cellDim;comp++) {
420 for (int x_order=0;x_order<basis_order;x_order++) {
421 for (int y_order=0;y_order<basis_order-x_order;y_order++) {
422 for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
423 fe_matrix.initialize();
424 // copy mass matrix
425 for (int i=0;i<numFields;i++) {
426 for (int j=0;j<numFields;j++) {
427 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
428 }
429 }
430
431 // clear old vector data
432 rhs_and_soln_vec.initialize();
433
434 // now get rhs vector
435 cub_points_cell.resize(1,numCubPointsCell,cellDim);
436
437 rhs_at_cub_points_cell.initialize();
438 rhsFunc(rhs_at_cub_points_cell,
439 cub_points_cell,
440 comp,
441 x_order,
442 y_order,
443 z_order);
444
445 cub_points_cell.resize(numCubPointsCell,cellDim);
446
447 cub_weights_cell.resize(numCubPointsCell);
448
449 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
450 rhs_at_cub_points_cell,
451 w_value_of_basis_at_cub_points_cell,
452 COMP_BLAS);
453
454 // now I need to get the boundary condition terms in place
455
456 for (unsigned i=0;i<4;i++) {
457 // map side quadrature to reference cell side
458 CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell,cub_points_side,sideDim,
459 (int) i, cell);
460
461 //evaluate basis at these points
462 basis->getValues( value_of_basis_at_cub_points_side_refcell ,
463 cub_points_side_refcell ,
464 OPERATOR_VALUE );
465
466 // evaluate imposed current on surface, which is n x curl(u_exact), on the quad points
467 cub_points_side_refcell.resize( 1 , numCubPointsSide , cellDim );
468 bcFunc(bc_func_at_cub_points_side_refcell,
469 cub_points_side_refcell,
470 cell ,
471 i,
472 comp ,
473 x_order,
474 y_order,
475 z_order);
476 cub_points_side_refcell.resize( numCubPointsSide , cellDim );
477
478 // now I need to integrate the bc function against the test basis
479 // need to weight the basis functions with quadrature weights
480 // the size of the face is embedded in the normal
481 cub_weights_side.resize(1,numCubPointsSide);
482 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_side_refcell ,
483 cub_weights_side ,
484 value_of_basis_at_cub_points_side_refcell );
485 cub_weights_side.resize(numCubPointsSide);
486
487 FunctionSpaceTools::integrate<double>( bc_fields_per_side ,
488 bc_func_at_cub_points_side_refcell ,
489 w_value_of_basis_at_cub_points_side_refcell ,
490 COMP_BLAS );
491
492 RealSpaceTools<double>::subtract(rhs_and_soln_vec, bc_fields_per_side );
493
494
495 }
496
497 // solve linear system
498 solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info);
499 solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
500
501 interp_points_ref.resize(1,numInterpPoints,cellDim);
502 // get exact solution for comparison
503 FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
504 exact_solution.initialize();
505 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
506 interp_points_ref.resize(numInterpPoints,cellDim);
507
508 // compute interpolant
509 // first evaluate basis at interpolation points
510 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
511 FunctionSpaceTools::evaluate<double>( interpolant ,
512 rhs_and_soln_vec ,
513 value_of_basis_at_interp_points );
514 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
515
516 RealSpaceTools<double>::subtract(interpolant,exact_solution);
517
518 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
519
520 *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
521 << x_order << ", " << y_order << ", " << z_order
522 << ") in component " << comp
523 << " and finite element interpolant of order " << basis_order << ": "
524 << nrm << "\n";
525
526 if (nrm > zero) {
527 *outStream << "\n\nPatch test failed for solution polynomial order ("
528 << x_order << ", " << y_order << ", " << z_order << ") and basis order "
529 << basis_order << "\n\n";
530 errorFlag++;
531 }
532 }
533 }
534 }
535 }
536 }
537
538 }
539
540 catch (const std::logic_error & err) {
541 *outStream << err.what() << "\n\n";
542 errorFlag = -1000;
543 };
544
545 if (errorFlag != 0)
546 std::cout << "End Result: TEST FAILED\n";
547 else
548 std::cout << "End Result: TEST PASSED\n";
549
550 // reset format state of std::cout
551 std::cout.copyfmt(oldFormatState);
552
553 return errorFlag;
554}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HCURL_TET_In_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getReferenceFaceNormal(ArrayFaceNormal &refFaceNormal, const int &faceOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to faces of 3D reference 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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
int dimension(const int whichDim) const
Returns the specified dimension.
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2.
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.