Intrepid
test_02.cpp
Go to the documentation of this file.
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
59#include "Teuchos_oblackholestream.hpp"
60#include "Teuchos_RCP.hpp"
61#include "Teuchos_GlobalMPISession.hpp"
62#include "Teuchos_SerialDenseMatrix.hpp"
63#include "Teuchos_SerialDenseVector.hpp"
64#include "Teuchos_LAPACK.hpp"
65
66using namespace std;
67using namespace Intrepid;
68
69void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int, int );
70void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int );
71
72void 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
88void rhsFunc( FieldContainer<double> & result ,
89 const FieldContainer<double> &points ,
90 int comp,
91 int xd,
92 int yd,
93 int zd )
94{
95 u_exact( result , points , comp , xd , yd , zd );
96}
97
98int main(int argc, char *argv[]) {
99 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
100
101 // This little trick lets us print to std::cout only if
102 // a (dummy) command-line argument is provided.
103 int iprint = argc - 1;
104 Teuchos::RCP<std::ostream> outStream;
105 Teuchos::oblackholestream bhs; // outputs nothing
106 if (iprint > 0)
107 outStream = Teuchos::rcp(&std::cout, false);
108 else
109 outStream = Teuchos::rcp(&bhs, false);
110
111 // Save the format state of the original std::cout.
112 Teuchos::oblackholestream oldFormatState;
113 oldFormatState.copyfmt(std::cout);
114
115 *outStream \
116 << "===============================================================================\n" \
117 << "| |\n" \
118 << "| Unit Test (Basis_HCURL_HEX_In_FEM) |\n" \
119 << "| |\n" \
120 << "| 1) Patch test involving H(curl) matrices |\n" \
121 << "| |\n" \
122 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
123 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
124 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
125 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
126 << "| |\n" \
127 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
128 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
129 << "| |\n" \
130 << "===============================================================================\n" \
131 << "| TEST 2: Patch test for mass matrices |\n" \
132 << "===============================================================================\n";
133
134
135 int errorFlag = 0;
136
137 outStream -> precision(16);
138
139 try {
140 shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<8> >()); // create parent cell topology
141
142 int cellDim = cell.getDimension();
143
144 int min_order = 1;
145 int max_order = 3;
146
147 int numIntervals = max_order;
148 int numInterpPoints = (numIntervals + 1)*(numIntervals +1)*(numIntervals+1);
149 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
150 int counter = 0;
151 for (int k=0; k<=numIntervals; k++) {
152 for (int j=0; j<=numIntervals; j++) {
153 for (int i=0;i<numIntervals;i++) {
154 interp_points_ref(counter,0) = i*(1.0/numIntervals);
155 interp_points_ref(counter,1) = j*(1.0/numIntervals);
156 interp_points_ref(counter,2) = k*(1.0/numIntervals);
157 counter++;
158 }
159 }
160 }
161
162 for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
163 // create basis
164 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
165 Teuchos::rcp(new Basis_HCURL_HEX_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_SPECTRAL) );
166 int numFields = basis->getCardinality();
167
168 // create cubatures
170 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create( cell , 2* basis_order );
171
172
173 int numCubPointsCell = cellCub->getNumPoints();
174
175 // hold cubature information
176 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
177 FieldContainer<double> cub_weights_cell(numCubPointsCell);
178
179 // hold basis function information on refcell
180 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
181 FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
182
183
184 // holds rhs data
185 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim);
186
187 // FEM mass matrix
188 FieldContainer<double> fe_matrix_bak(1,numFields,numFields);
189 FieldContainer<double> fe_matrix(1,numFields,numFields);
190 FieldContainer<double> rhs_and_soln_vec(1,numFields);
191
192 FieldContainer<int> ipiv(numFields);
193 FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim);
194 FieldContainer<double> interpolant( 1, numInterpPoints , cellDim );
195
196 int info = 0;
197 Teuchos::LAPACK<int, double> solver;
198
199 // set test tolerance
200 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
201
202 // build matrices outside the loop, and then just do the rhs
203 // for each iteration
204 cellCub->getCubature(cub_points_cell, cub_weights_cell);
205
206 // need the vector basis
207 basis->getValues(value_of_basis_at_cub_points_cell,
208 cub_points_cell,
209 OPERATOR_VALUE);
210 basis->getValues( value_of_basis_at_interp_points ,
211 interp_points_ref ,
212 OPERATOR_VALUE );
213
214
215 // construct mass matrix
216 cub_weights_cell.resize(1,numCubPointsCell);
217 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
218 cub_weights_cell ,
219 value_of_basis_at_cub_points_cell );
220 cub_weights_cell.resize(numCubPointsCell);
221
222
223 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
224 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
225 w_value_of_basis_at_cub_points_cell ,
226 value_of_basis_at_cub_points_cell ,
227 COMP_BLAS );
228 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
229
230 for (int x_order=0;x_order<basis_order;x_order++) {
231 for (int y_order=0;y_order<basis_order;y_order++) {
232 for (int z_order=0;z_order<basis_order;z_order++) {
233 for (int comp=0;comp<cellDim;comp++) {
234 fe_matrix.initialize();
235 // copy mass matrix
236 for (int i=0;i<numFields;i++) {
237 for (int j=0;j<numFields;j++) {
238 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
239 }
240 }
241
242 // clear old vector data
243 rhs_and_soln_vec.initialize();
244
245 // now get rhs vector
246
247 cub_points_cell.resize(1,numCubPointsCell,cellDim);
248
249 rhs_at_cub_points_cell.initialize();
250 rhsFunc(rhs_at_cub_points_cell,
251 cub_points_cell,
252 comp,
253 x_order,
254 y_order,
255 z_order);
256
257 cub_points_cell.resize(numCubPointsCell,cellDim);
258 cub_weights_cell.resize(numCubPointsCell);
259
260 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
261 rhs_at_cub_points_cell,
262 w_value_of_basis_at_cub_points_cell,
263 COMP_BLAS);
264
265 // solve linear system
266
267 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vec[0],
268 numFields, &info);
269// solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info);
270// solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
271
272 interp_points_ref.resize(1,numInterpPoints,cellDim);
273 // get exact solution for comparison
274 FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
275 exact_solution.initialize();
276 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
277 interp_points_ref.resize(numInterpPoints,cellDim);
278
279 // compute interpolant
280 // first evaluate basis at interpolation points
281 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
282 FunctionSpaceTools::evaluate<double>( interpolant ,
283 rhs_and_soln_vec ,
284 value_of_basis_at_interp_points );
285 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
286
287 RealSpaceTools<double>::subtract(interpolant,exact_solution);
288
289 // let's compute the L2 norm
290 interpolant.resize(1,numInterpPoints,cellDim);
291 FieldContainer<double> l2norm(1);
292 FunctionSpaceTools::dataIntegral<double>( l2norm ,
293 interpolant ,
294 interpolant ,
295 COMP_BLAS );
296
297 const double nrm = sqrtl(l2norm(0));
298
299 *outStream << "\nFunction space L^2 Norm-2 error between scalar components of exact solution of order ("
300 << x_order << ", " << y_order << ", " << z_order
301 << ") in component " << comp
302 << " and finite element interpolant of order " << basis_order << ": "
303 << nrm << "\n";
304
305 if (nrm > zero) {
306 *outStream << "\n\nPatch test failed for solution polynomial order ("
307 << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector) ("
308 << basis_order << ", " << basis_order+1 << ")\n\n";
309 errorFlag++;
310 }
311 }
312 }
313 }
314 }
315 }
316
317 }
318
319 catch (const std::logic_error & err) {
320 *outStream << err.what() << "\n\n";
321 errorFlag = -1000;
322 };
323
324 if (errorFlag != 0)
325 std::cout << "End Result: TEST FAILED\n";
326 else
327 std::cout << "End Result: TEST PASSED\n";
328
329 // reset format state of std::cout
330 std::cout.copyfmt(oldFormatState);
331
332 return errorFlag;
333}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::CubatureTensor class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HCURL_HEX_In_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell.
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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int dimension(const int whichDim) const
Returns the specified dimension.
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2.