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 
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, int );
67 void bcFunc( FieldContainer<double> &, const FieldContainer<double> & ,
68  const shards::CellTopology & ,
69  int , int , int , int , int );
70 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int );
71 
72 void 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 
88 void 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 
170 void 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 
259 int 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 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for the Intrepid::HCURL_TET_In_FEM class.
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 the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
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: