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