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