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