Intrepid
test_04.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_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 #define INTREPID_TEST_COMMAND( S ) \
63 { \
64  try { \
65  S ; \
66  } \
67  catch (std::logic_error err) { \
68  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
69  *outStream << err.what() << '\n'; \
70  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
71  }; \
72 }
73 
74 
75 int main(int argc, char *argv[]) {
76 
77  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
78 
79  // This little trick lets us print to std::cout only if
80  // a (dummy) command-line argument is provided.
81  int iprint = argc - 1;
82  Teuchos::RCP<std::ostream> outStream;
83  Teuchos::oblackholestream bhs; // outputs nothing
84  if (iprint > 0)
85  outStream = Teuchos::rcp(&std::cout, false);
86  else
87  outStream = Teuchos::rcp(&bhs, false);
88 
89  // Save the format state of the original std::cout.
90  Teuchos::oblackholestream oldFormatState;
91  oldFormatState.copyfmt(std::cout);
92 
93  *outStream \
94  << "===============================================================================\n" \
95  << "| |\n" \
96  << "| Unit Test (FunctionSpaceTools) |\n" \
97  << "| |\n" \
98  << "| 1) volume integration on tetrahedra, testing dataIntegral |\n" \
99  << "| |\n" \
100  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
101  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
102  << "| |\n" \
103  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105  << "| |\n" \
106  << "===============================================================================\n";
107 
108 
109  int errorFlag = 0;
110 
111  typedef FunctionSpaceTools fst;
112 
113  *outStream \
114  << "\n"
115  << "===============================================================================\n"\
116  << "| TEST 1: correctness of cell volumes |\n"\
117  << "===============================================================================\n";
118 
119  outStream->precision(20);
120 
121  try {
122  shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tet
123 
124  /* Related to cubature. */
125  DefaultCubatureFactory<double> cubFactory; // create cubature factory
126  int cubDegree = 0; // cubature degree
127  Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature
128  int spaceDim = myCub->getDimension(); // get spatial dimension
129  int numCubPoints = myCub->getNumPoints(); // get number of cubature points
130 
131  /* Cell geometries. */
132  int numCells = 4;
133  int numNodes = 4;
134  int numCellData = numCells*numNodes*spaceDim;
135  double tetnodes[] = {
136  // tet 0
137  0.0, 0.0, 0.0,
138  1.0, 0.0, 0.0,
139  0.0, 1.0, 0.0,
140  0.0, 0.0, 1.0,
141  // tet 1
142  4.0, 5.0, 1.0,
143  -6.0, 2.0, 0.0,
144  4.0, -3.0, -1.0,
145  0.0, 2.0, 5.0,
146  // tet 2
147  -6.0, -3.0, 1.0,
148  9.0, 2.0, 1.0,
149  8.9, 2.1, 0.9,
150  8.9, 2.1, 1.1,
151  // tet 3
152  -6.0, -3.0, 1.0,
153  12.0, 3.0, 1.0,
154  2.9, 0.1, 0.9,
155  2.9, 0.1, 1.1
156  };
157 
158  /* Analytic volumes. */
159  double tetvols[] = {1.0/6.0, 194.0/3.0, 1.0/15.0, 2.0/25.0};
160 
161  /* Computational arrays. */
162  FieldContainer<double> cub_points(numCubPoints, spaceDim);
163  FieldContainer<double> cub_weights(numCubPoints);
164  FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
165  FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
166  FieldContainer<double> jacobian_det(numCells, numCubPoints);
167  FieldContainer<double> weighted_measure(numCells, numCubPoints);
168  FieldContainer<double> data_one(numCells, numCubPoints);
169  FieldContainer<double> volumes(numCells);
170 
171  /******************* START COMPUTATION ***********************/
172 
173  // get cubature points and weights
174  myCub->getCubature(cub_points, cub_weights);
175 
176  // fill cell vertex array
177  cell_nodes.setValues(tetnodes, numCellData);
178 
179  // compute geometric cell information
180  CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
181  CellTools<double>::setJacobianDet(jacobian_det, jacobian);
182 
183  // compute weighted measure
184  fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
185 
186  // set data to 1.0
187  for (int cell=0; cell<data_one.dimension(0); cell++) {
188  for (int qp=0; qp<data_one.dimension(1); qp++) {
189  data_one(cell,qp) = 1.0;
190  }
191  }
192 
193  // compute volumes
194  fst::integrate<double>(volumes, data_one, weighted_measure, COMP_CPP);
195 
196  /******************* STOP COMPUTATION ***********************/
197 
198 
199  /******************* START COMPARISON ***********************/
200  for (int cell_id = 0; cell_id < numCells; cell_id++) {
201  *outStream << "Volume of cell " << cell_id << " = " << volumes(cell_id) << " vs. Analytic value = " << tetvols[cell_id] << "\n";
202  if (std::fabs(volumes(cell_id)-tetvols[cell_id]) > INTREPID_TOL) {
203  errorFlag++;
204  }
205  }
206  /******************* STOP COMPARISON ***********************/
207 
208  *outStream << "\n";
209  }
210  catch (std::logic_error err) {
211  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
212  *outStream << err.what() << '\n';
213  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
214  errorFlag = -1000;
215  };
216 
217 
218  if (errorFlag != 0)
219  std::cout << "End Result: TEST FAILED\n";
220  else
221  std::cout << "End Result: TEST PASSED\n";
222 
223  // reset format state of std::cout
224  std::cout.copyfmt(oldFormatState);
225 
226  return errorFlag;
227 }
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...
Intrepid utilities.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
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: