Intrepid
test_05.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 
52 #include "Intrepid_FunctionSpaceTools.hpp"
53 #include "Intrepid_FieldContainer.hpp"
54 #include "Shards_Array.hpp"
55 #include "Intrepid_CellTools.hpp"
56 #include "Intrepid_HDIV_HEX_I1_FEM.hpp"
57 #include "Intrepid_DefaultCubatureFactory.hpp"
58 #include "Intrepid_Utils.hpp"
59 #include "Teuchos_oblackholestream.hpp"
60 #include "Teuchos_RCP.hpp"
61 #include "Teuchos_GlobalMPISession.hpp"
62 
63 using namespace std;
64 using namespace Intrepid;
65 
66 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell )
67 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell )
68 
69 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field )
70 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field )
71 
72 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point )
73 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point )
74 
75 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim )
76 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim )
77 
78 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Node )
79 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Node )
80 
81 
82 typedef FieldContainer<double> FC;
83 
84 
85 int main(int argc, char *argv[]) {
86 
87  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
88 
89  // This little trick lets us print to std::cout only if
90  // a (dummy) command-line argument is provided.
91  int iprint = argc - 1;
92  Teuchos::RCP<std::ostream> outStream;
93  Teuchos::oblackholestream bhs; // outputs nothing
94  if (iprint > 0)
95  outStream = Teuchos::rcp(&std::cout, false);
96  else
97  outStream = Teuchos::rcp(&bhs, false);
98 
99  // Save the format state of the original std::cout.
100  Teuchos::oblackholestream oldFormatState;
101  oldFormatState.copyfmt(std::cout);
102 
103  *outStream \
104  << "===============================================================================\n" \
105  << "| |\n" \
106  << "| Unit Test (FunctionSpaceTools) |\n" \
107  << "| |\n" \
108  << "| 1) basic operator transformations and integration in HDIV |\n" \
109  << "| via shards::Array wrappers |\n" \
110  << "| |\n" \
111  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
112  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
113  << "| |\n" \
114  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
115  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
116  << "| |\n" \
117  << "===============================================================================\n";
118 
119 
120  int errorFlag = 0;
121 
122  typedef FunctionSpaceTools fst;
123 
124  *outStream \
125  << "\n"
126  << "===============================================================================\n"\
127  << "| TEST 1: correctness of math operations |\n"\
128  << "===============================================================================\n";
129 
130  outStream->precision(20);
131 
132  try {
133  shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >(); // cell type: hex
134 
135  /* Related to cubature. */
136  DefaultCubatureFactory<double> cubFactory; // create cubature factory
137  int cubDegree = 20; // cubature degree
138  Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature
139  int spaceDim = myCub->getDimension(); // get spatial dimension
140  int numCubPoints = myCub->getNumPoints(); // get number of cubature points
141 
142  /* Related to basis. */
143  Basis_HDIV_HEX_I1_FEM<double, FieldContainer<double> > hexBasis; // create H-div basis on a hex
144  int numFields = hexBasis.getCardinality(); // get basis cardinality
145 
146  /* Cell geometries and orientations. */
147  int numCells = 4;
148  int numNodes = 8;
149  int numCellData = numCells*numNodes*spaceDim;
150  int numSignData = numCells*numFields;
151 
152  double hexnodes[] = {
153  // hex 0 -- affine
154  -1.0, -1.0, -1.0,
155  1.0, -1.0, -1.0,
156  1.0, 1.0, -1.0,
157  -1.0, 1.0, -1.0,
158  -1.0, -1.0, 1.0,
159  1.0, -1.0, 1.0,
160  1.0, 1.0, 1.0,
161  -1.0, 1.0, 1.0,
162  // hex 1 -- affine
163  -3.0, -3.0, 1.0,
164  6.0, 3.0, 1.0,
165  7.0, 8.0, 0.0,
166  -2.0, 2.0, 0.0,
167  -3.0, -3.0, 4.0,
168  6.0, 3.0, 4.0,
169  7.0, 8.0, 3.0,
170  -2.0, 2.0, 3.0,
171  // hex 2 -- affine
172  -3.0, -3.0, 0.0,
173  9.0, 3.0, 0.0,
174  15.0, 6.1, 0.0,
175  3.0, 0.1, 0.0,
176  9.0, 3.0, 0.1,
177  21.0, 9.0, 0.1,
178  27.0, 12.1, 0.1,
179  15.0, 6.1, 0.1,
180  // hex 3 -- nonaffine
181  -2.0, -2.0, 0.0,
182  2.0, -1.0, 0.0,
183  1.0, 6.0, 0.0,
184  -1.0, 1.0, 0.0,
185  0.0, 0.0, 1.0,
186  1.0, 0.0, 1.0,
187  1.0, 1.0, 1.0,
188  0.0, 1.0, 1.0
189  };
190 
191  short facesigns[] = {
192  1, 1, 1, 1, 1, 1,
193  1, -1, 1, -1, 1, -1,
194  -1, -1, 1, 1, -1, 1,
195  -1, -1, 1, 1, -1, -1
196  };
197 
198  /* Computational arrays. */
199  // First allocate one very large work space.
200  Teuchos::Array<double> work_space(numCubPoints*spaceDim +
201  numCubPoints +
202  numCells*numNodes*spaceDim +
203  numCells*numCubPoints*spaceDim*spaceDim +
204  2*numCells*numCubPoints +
205  numFields*numCubPoints +
206  2*numCells*numFields*numCubPoints +
207  numCells*numFields*numFields +
208  numFields*numCubPoints*spaceDim +
209  2*numCells*numFields*numCubPoints*spaceDim +
210  numCells*numFields*numFields
211  );
212 
213  int offset = 0;
214  shards::Array<double,shards::NaturalOrder,Point,Dim> cub_points(&work_space[offset], numCubPoints, spaceDim);
215  FC cub_points_FC(cub_points);
216  offset += numCubPoints*spaceDim;
217  shards::Array<double,shards::NaturalOrder,Point> cub_weights(&work_space[offset], numCubPoints);
218  FC cub_weights_FC(cub_weights);
219  offset += numCubPoints;
220  shards::Array<double,shards::NaturalOrder,Cell,Node,Dim> cell_nodes(&work_space[offset], numCells, numNodes, spaceDim);
221  FC cell_nodes_FC(cell_nodes);
222  offset += numCells*numNodes*spaceDim;
223  FieldContainer<short> field_signs(numCells, numFields);
224  shards::Array<double,shards::NaturalOrder,Cell,Point,Dim,Dim> jacobian(&work_space[offset], numCells, numCubPoints, spaceDim, spaceDim);
225  FC jacobian_FC(jacobian);
226  offset += numCells*numCubPoints*spaceDim*spaceDim;
227  //shards::Array<double,shards::NaturalOrder> jacobian_inv(&work_space[offset]numCells, numCubPoints, spaceDim, spaceDim);
228  shards::Array<double,shards::NaturalOrder,Cell,Point> jacobian_det(&work_space[offset], numCells, numCubPoints);
229  FC jacobian_det_FC(jacobian_det);
230  offset += numCells*numCubPoints;
231  shards::Array<double,shards::NaturalOrder,Cell,Point> weighted_measure(&work_space[offset], numCells, numCubPoints);
232  FC weighted_measure_FC(weighted_measure);
233  offset += numCells*numCubPoints;
234 
235  shards::Array<double,shards::NaturalOrder,Field,Point> div_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints);
236  FC div_of_basis_at_cub_points_FC(div_of_basis_at_cub_points);
237  offset += numFields*numCubPoints;
238  shards::Array<double,shards::NaturalOrder,Cell,Field,Point>
239  transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints);
240  FC transformed_div_of_basis_at_cub_points_FC(transformed_div_of_basis_at_cub_points);
241  offset += numCells*numFields*numCubPoints;
242  shards::Array<double,shards::NaturalOrder,Cell,Field,Point>
243  weighted_transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints);
244  FC weighted_transformed_div_of_basis_at_cub_points_FC(weighted_transformed_div_of_basis_at_cub_points);
245  offset += numCells*numFields*numCubPoints;
246  shards::Array<double,shards::NaturalOrder,Cell,Field,Field> stiffness_matrices(&work_space[offset], numCells, numFields, numFields);
247  FC stiffness_matrices_FC(stiffness_matrices);
248  offset += numCells*numFields*numFields;
249 
250  shards::Array<double,shards::NaturalOrder,Field,Point,Dim> value_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints, spaceDim);
251  FC value_of_basis_at_cub_points_FC(value_of_basis_at_cub_points);
252  offset += numFields*numCubPoints*spaceDim;
253  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim>
254  transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim);
255  FC transformed_value_of_basis_at_cub_points_FC(transformed_value_of_basis_at_cub_points);
256  offset += numCells*numFields*numCubPoints*spaceDim;
257  shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim>
258  weighted_transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim);
259  FC weighted_transformed_value_of_basis_at_cub_points_FC(weighted_transformed_value_of_basis_at_cub_points);
260  offset += numCells*numFields*numCubPoints*spaceDim;
261  shards::Array<double,shards::NaturalOrder,Cell,Field,Field> mass_matrices(&work_space[offset], numCells, numFields, numFields);
262  FC mass_matrices_FC(mass_matrices);
263 
264 
265  /******************* START COMPUTATION ***********************/
266 
267  // get cubature points and weights
268  myCub->getCubature(cub_points_FC, cub_weights_FC);
269 
270  // fill cell vertex array
271  cell_nodes_FC.setValues(hexnodes, numCellData);
272 
273  // set basis function signs, for each cell
274  field_signs.setValues(facesigns, numSignData);
275 
276  // compute geometric cell information
277  CellTools<double>::setJacobian(jacobian_FC, cub_points_FC, cell_nodes_FC, cellType);
278  //CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
279  CellTools<double>::setJacobianDet(jacobian_det_FC, jacobian_FC);
280 
281  // compute weighted measure
282  fst::computeCellMeasure<double>(weighted_measure_FC, jacobian_det_FC, cub_weights_FC);
283 
284 
285  // Computing stiffness matrices:
286  // tabulate divergences of basis functions at (reference) cubature points
287  hexBasis.getValues(div_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_DIV);
288 
289  // transform divergences of basis functions
290  fst::HDIVtransformDIV<double>(transformed_div_of_basis_at_cub_points_FC,
291  jacobian_det_FC,
292  div_of_basis_at_cub_points_FC);
293 
294  // multiply with weighted measure
295  fst::multiplyMeasure<double>(weighted_transformed_div_of_basis_at_cub_points_FC,
296  weighted_measure_FC,
297  transformed_div_of_basis_at_cub_points_FC);
298 
299  // compute stiffness matrices
300  fst::integrate<double>(stiffness_matrices_FC,
301  transformed_div_of_basis_at_cub_points_FC,
302  weighted_transformed_div_of_basis_at_cub_points_FC,
303  COMP_CPP);
304 
305  // apply field signs
306  fst::applyLeftFieldSigns<double>(stiffness_matrices_FC, field_signs);
307  fst::applyRightFieldSigns<double>(stiffness_matrices_FC, field_signs);
308 
309 
310  // Computing mass matrices:
311  // tabulate values of basis functions at (reference) cubature points
312  hexBasis.getValues(value_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_VALUE);
313 
314  // transform values of basis functions
315  fst::HDIVtransformVALUE<double>(transformed_value_of_basis_at_cub_points_FC,
316  jacobian_FC,
317  jacobian_det_FC,
318  value_of_basis_at_cub_points_FC);
319 
320  // multiply with weighted measure
321  fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_FC,
322  weighted_measure_FC,
323  transformed_value_of_basis_at_cub_points_FC);
324 
325  // compute mass matrices
326  fst::integrate<double>(mass_matrices_FC,
327  transformed_value_of_basis_at_cub_points_FC,
328  weighted_transformed_value_of_basis_at_cub_points_FC,
329  COMP_CPP);
330 
331  // apply field signs
332  fst::applyLeftFieldSigns<double>(mass_matrices_FC, field_signs);
333  fst::applyRightFieldSigns<double>(mass_matrices_FC, field_signs);
334 
335  /******************* STOP COMPUTATION ***********************/
336 
337 
338  /******************* START COMPARISON ***********************/
339  string basedir = "./testdata";
340  for (int cell_id = 0; cell_id < numCells-1; cell_id++) {
341 
342  stringstream namestream;
343  string filename;
344  namestream << basedir << "/mass_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
345  namestream >> filename;
346 
347  ifstream massfile(&filename[0]);
348  if (massfile.is_open()) {
349  if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
350  errorFlag++;
351  massfile.close();
352  }
353  else {
354  errorFlag = -1;
355  std::cout << "End Result: TEST FAILED\n";
356  return errorFlag;
357  }
358 
359  namestream.clear();
360  namestream << basedir << "/stiff_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
361  namestream >> filename;
362 
363  ifstream stifffile(&filename[0]);
364  if (stifffile.is_open())
365  {
366  if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
367  errorFlag++;
368  stifffile.close();
369  }
370  else {
371  errorFlag = -1;
372  std::cout << "End Result: TEST FAILED\n";
373  return errorFlag;
374  }
375 
376  }
377  for (int cell_id = 3; cell_id < numCells; cell_id++) {
378 
379  stringstream namestream;
380  string filename;
381  namestream << basedir << "/mass_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
382  namestream >> filename;
383 
384  ifstream massfile(&filename[0]);
385  if (massfile.is_open()) {
386  if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
387  errorFlag++;
388  massfile.close();
389  }
390  else {
391  errorFlag = -1;
392  std::cout << "End Result: TEST FAILED\n";
393  return errorFlag;
394  }
395 
396  namestream.clear();
397  namestream << basedir << "/stiff_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
398  namestream >> filename;
399 
400  ifstream stifffile(&filename[0]);
401  if (stifffile.is_open())
402  {
403  if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
404  errorFlag++;
405  stifffile.close();
406  }
407  else {
408  errorFlag = -1;
409  std::cout << "End Result: TEST FAILED\n";
410  return errorFlag;
411  }
412 
413  }
414  /******************* STOP COMPARISON ***********************/
415 
416  *outStream << "\n";
417  }
418  catch (std::logic_error err) {
419  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
420  *outStream << err.what() << '\n';
421  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
422  errorFlag = -1000;
423  };
424 
425 
426  if (errorFlag != 0)
427  std::cout << "End Result: TEST FAILED\n";
428  else
429  std::cout << "End Result: TEST PASSED\n";
430 
431  // reset format state of std::cout
432  std::cout.copyfmt(oldFormatState);
433 
434  return errorFlag;
435 }
Intrepid::Basis_HDIV_HEX_I1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HDIV_HEX_I1_FEMDef.hpp:95
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::FunctionSpaceTools
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
Definition: Intrepid_FunctionSpaceTools.hpp:69
Intrepid_HDIV_HEX_I1_FEM.hpp
Header file for the Intrepid::HDIV_HEX_I1_FEM class.
Intrepid_Utils.hpp
Intrepid utilities.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::Basis_HDIV_HEX_I1_FEM
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedron cell...
Definition: Intrepid_HDIV_HEX_I1_FEM.hpp:105
Intrepid::FieldContainer
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Definition: Intrepid_FieldContainer.hpp:78
Intrepid_FunctionSpaceTools.hpp
Header file for the Intrepid::FunctionSpaceTools class.
Intrepid::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition: Intrepid_DefaultCubatureFactory.hpp:77
Intrepid::DefaultCubatureFactory::create
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Definition: Intrepid_DefaultCubatureFactoryDef.hpp:53
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:111
Generated by   doxygen 1.8.5