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
58#include "Teuchos_oblackholestream.hpp"
59#include "Teuchos_RCP.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
61#include "Teuchos_SerialDenseMatrix.hpp"
62#include "Teuchos_SerialDenseVector.hpp"
63#include "Teuchos_LAPACK.hpp"
64
65using namespace std;
66using namespace Intrepid;
67
68void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int );
69void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int, int );
70
71// This is the rhs for (div tau,w) = (f,w),
72// which makes f the negative Laplacian of scalar solution
74 const FieldContainer<double> &points,
75 int xd,
76 int yd )
77{
78 for (int cell=0;cell<result.dimension(0);cell++) {
79 for (int pt=0;pt<result.dimension(1);pt++) {
80 result(cell,pt) = 0.0;
81 if (xd >=2) {
82 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd);
83 }
84 if (yd >=2) {
85 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2);
86 }
87 }
88 }
89}
90
92 const FieldContainer<double> &points,
93 int xd,
94 int yd)
95{
96 for (int cell=0;cell<result.dimension(0);cell++){
97 for (int pt=0;pt<result.dimension(1);pt++) {
98 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd);
99 }
100 }
101 return;
102}
103
104int main(int argc, char *argv[]) {
105 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
106
107 // This little trick lets us print to std::cout only if
108 // a (dummy) command-line argument is provided.
109 int iprint = argc - 1;
110 Teuchos::RCP<std::ostream> outStream;
111 Teuchos::oblackholestream bhs; // outputs nothing
112 if (iprint > 0)
113 outStream = Teuchos::rcp(&std::cout, false);
114 else
115 outStream = Teuchos::rcp(&bhs, false);
116
117 // Save the format state of the original std::cout.
118 Teuchos::oblackholestream oldFormatState;
119 oldFormatState.copyfmt(std::cout);
120
121 *outStream \
122 << "===============================================================================\n" \
123 << "| |\n" \
124 << "| Unit Test (Basis_HGRAD_QUAD_Cn_FEM) |\n" \
125 << "| |\n" \
126 << "| 1) Patch test involving H(div) matrices |\n" \
127 << "| for the Dirichlet problem on a triangular patch |\n" \
128 << "| Omega with boundary Gamma. |\n" \
129 << "| |\n" \
130 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
131 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
132 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
133 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
134 << "| |\n" \
135 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
136 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
137 << "| |\n" \
138 << "===============================================================================\n"\
139 << "| TEST 1: Patch test |\n"\
140 << "===============================================================================\n";
141
142
143 int errorFlag = 0;
144
145 outStream -> precision(16);
146
147 try {
149 shards::CellTopology cell(shards::getCellTopologyData< shards::Quadrilateral<> >());
150 shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >());
151
152 int cellDim = cell.getDimension();
153 int sideDim = side.getDimension();
154
155 int min_order = 0;
156 int max_order = 5;
157
158 int numIntervals = max_order;
159 int numInterpPoints = (numIntervals + 1)*(numIntervals + 1);;
160 FieldContainer<double> interp_points_ref(numInterpPoints, 2);
161 int counter = 0;
162 for (int j=0; j<=numIntervals; j++) {
163 for (int i=0; i<=numIntervals; i++) {
164 interp_points_ref(counter,0) = i*(1.0/numIntervals);
165 interp_points_ref(counter,1) = j*(1.0/numIntervals);
166 counter++;
167 }
168 }
169
170 interp_points_ref.resize(numInterpPoints,2);
171
172 for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
173 // create bases
174 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
175 Teuchos::rcp(new Basis_HDIV_QUAD_In_FEM<double,FieldContainer<double> >(basis_order+1, POINTTYPE_SPECTRAL) );
176
177 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
178 Teuchos::rcp(new Basis_HGRAD_QUAD_Cn_FEM<double,FieldContainer<double> >(basis_order, POINTTYPE_SPECTRAL) );
179
180 int numVectorFields = vectorBasis->getCardinality();
181 int numScalarFields = scalarBasis->getCardinality();
182 int numTotalFields = numVectorFields + numScalarFields;
183
184 // create cubatures
185 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
186 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
187
188 int numCubPointsCell = cellCub->getNumPoints();
189 int numCubPointsSide = sideCub->getNumPoints();
190
191 // hold cubature information
192 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
193 FieldContainer<double> cub_weights_cell(numCubPointsCell);
194 FieldContainer<double> cub_points_side( numCubPointsSide, sideDim );
195 FieldContainer<double> cub_weights_side( numCubPointsSide );
196 FieldContainer<double> cub_points_side_refcell( numCubPointsSide , cellDim );
197
198 // hold basis function information on refcell
199 FieldContainer<double> value_of_v_basis_at_cub_points_cell(numVectorFields, numCubPointsCell, cellDim );
200 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
201 FieldContainer<double> div_of_v_basis_at_cub_points_cell( numVectorFields, numCubPointsCell );
202 FieldContainer<double> w_div_of_v_basis_at_cub_points_cell( 1, numVectorFields , numCubPointsCell );
203 FieldContainer<double> value_of_s_basis_at_cub_points_cell(numScalarFields,numCubPointsCell);
204 FieldContainer<double> w_value_of_s_basis_at_cub_points_cell(1,numScalarFields,numCubPointsCell);
205
206 // containers for side integration:
207 // I just need the normal component of the vector basis
208 // and the exact solution at the cub points
209 FieldContainer<double> value_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide,cellDim);
210 FieldContainer<double> n_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide);
211 FieldContainer<double> w_n_of_v_basis_at_cub_points_side(1,numVectorFields,numCubPointsSide);
212 FieldContainer<double> diri_data_at_cub_points_side(1,numCubPointsSide);
213 FieldContainer<double> side_normal(cellDim);
214
215 // holds rhs data
216 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell);
217
218 // FEM matrices and vectors
219 FieldContainer<double> fe_matrix_M(1,numVectorFields,numVectorFields);
220 FieldContainer<double> fe_matrix_B(1,numVectorFields,numScalarFields);
221 FieldContainer<double> fe_matrix(1,numTotalFields,numTotalFields);
222
223 FieldContainer<double> rhs_vector_vec(1,numVectorFields);
224 FieldContainer<double> rhs_vector_scal(1,numScalarFields);
225 FieldContainer<double> rhs_and_soln_vec(1,numTotalFields);
226
227 FieldContainer<int> ipiv(numTotalFields);
228 FieldContainer<double> value_of_s_basis_at_interp_points( numScalarFields , numInterpPoints);
229 FieldContainer<double> interpolant( 1 , numInterpPoints );
230
231 // set test tolerance
232 double zero = (basis_order+1)*(basis_order+1)*100*INTREPID_TOL;
233
234 // build matrices outside the loop, and then just do the rhs
235 // for each iteration
236
237 cellCub->getCubature(cub_points_cell, cub_weights_cell);
238 sideCub->getCubature(cub_points_side, cub_weights_side);
239
240 // need the vector basis & its divergences
241 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
242 cub_points_cell,
243 OPERATOR_VALUE);
244 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
245 cub_points_cell,
246 OPERATOR_DIV);
247
248 // need the scalar basis as well
249 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
250 cub_points_cell,
251 OPERATOR_VALUE);
252
253 // construct mass matrix
254 cub_weights_cell.resize(1,numCubPointsCell);
255 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
256 cub_weights_cell ,
257 value_of_v_basis_at_cub_points_cell );
258 cub_weights_cell.resize(numCubPointsCell);
259
260
261 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
262 FunctionSpaceTools::integrate<double>(fe_matrix_M,
263 w_value_of_v_basis_at_cub_points_cell ,
264 value_of_v_basis_at_cub_points_cell ,
265 COMP_BLAS );
266 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
267
268 // div matrix
269 cub_weights_cell.resize(1,numCubPointsCell);
270 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
271 cub_weights_cell,
272 div_of_v_basis_at_cub_points_cell);
273 cub_weights_cell.resize(numCubPointsCell);
274
275 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
276 FunctionSpaceTools::integrate<double>(fe_matrix_B,
277 w_div_of_v_basis_at_cub_points_cell ,
278 value_of_s_basis_at_cub_points_cell ,
279 COMP_BLAS );
280 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
281
282
283 // construct div matrix
284
285 for (int x_order=0;x_order<=basis_order;x_order++) {
286 for (int y_order=0;y_order<=basis_order;y_order++) {
287
288 // reset global matrix since I destroyed it in LU factorization.
289 fe_matrix.initialize();
290 // insert mass matrix into global matrix
291 for (int i=0;i<numVectorFields;i++) {
292 for (int j=0;j<numVectorFields;j++) {
293 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
294 }
295 }
296
297 // insert div matrix into global matrix
298 for (int i=0;i<numVectorFields;i++) {
299 for (int j=0;j<numScalarFields;j++) {
300 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
301 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
302 }
303 }
304
305 // clear old vector data
306 rhs_vector_vec.initialize();
307 rhs_vector_scal.initialize();
308 rhs_and_soln_vec.initialize();
309
310 // now get rhs vector
311 // rhs_vector_scal is just (rhs,w) for w in the scalar basis
312 // I already have the scalar basis tabulated.
313 cub_points_cell.resize(1,numCubPointsCell,cellDim);
314 rhsFunc(rhs_at_cub_points_cell,
315 cub_points_cell,
316 x_order,
317 y_order);
318
319 cub_points_cell.resize(numCubPointsCell,cellDim);
320
321 cub_weights_cell.resize(1,numCubPointsCell);
322 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
323 cub_weights_cell,
324 value_of_s_basis_at_cub_points_cell);
325 cub_weights_cell.resize(numCubPointsCell);
326 FunctionSpaceTools::integrate<double>(rhs_vector_scal,
327 rhs_at_cub_points_cell,
328 w_value_of_s_basis_at_cub_points_cell,
329 COMP_BLAS);
330
331 for (int i=0;i<numScalarFields;i++) {
332 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
333 }
334
335
336 // now get <u,v.n> on boundary
337 for (unsigned side_cur=0;side_cur<4;side_cur++) {
338 // map side cubature to current side
339 CellTools<double>::mapToReferenceSubcell( cub_points_side_refcell ,
340 cub_points_side ,
341 sideDim ,
342 (int)side_cur ,
343 cell );
344
345 // Evaluate dirichlet data
346 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
347 u_exact(diri_data_at_cub_points_side,
348 cub_points_side_refcell,x_order,y_order);
349 cub_points_side_refcell.resize(numCubPointsSide,cellDim);
350
351 // get normal direction, this has the edge weight factored into it already
353 (int)side_cur,cell );
354
355 // v.n at cub points on side
356 vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
357 cub_points_side_refcell ,
358 OPERATOR_VALUE );
359
360
361 for (int i=0;i<numVectorFields;i++) {
362 for (int j=0;j<numCubPointsSide;j++) {
363 n_of_v_basis_at_cub_points_side(i,j) = 0.0;
364 for (int k=0;k<cellDim;k++) {
365 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) *
366 value_of_v_basis_at_cub_points_side(i,j,k);
367 }
368 }
369 }
370
371 cub_weights_side.resize(1,numCubPointsSide);
372 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
373 cub_weights_side,
374 n_of_v_basis_at_cub_points_side);
375 cub_weights_side.resize(numCubPointsSide);
376
377 FunctionSpaceTools::integrate<double>(rhs_vector_vec,
378 diri_data_at_cub_points_side,
379 w_n_of_v_basis_at_cub_points_side,
380 COMP_BLAS,
381 false);
382 for (int i=0;i<numVectorFields;i++) {
383 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
384 }
385
386 }
387
388 // solve linear system
389 int info = 0;
390 Teuchos::LAPACK<int, double> solver;
391 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0],
392 numTotalFields, &info);
393
394 // compute interpolant; the scalar entries are last
395 scalarBasis->getValues(value_of_s_basis_at_interp_points,
396 interp_points_ref,
397 OPERATOR_VALUE);
398 for (int pt=0;pt<numInterpPoints;pt++) {
399 interpolant(0,pt)=0.0;
400 for (int i=0;i<numScalarFields;i++) {
401 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
402 * value_of_s_basis_at_interp_points(i,pt);
403 }
404 }
405
406 interp_points_ref.resize(1,numInterpPoints,cellDim);
407 // get exact solution for comparison
408 FieldContainer<double> exact_solution(1,numInterpPoints);
409 u_exact( exact_solution , interp_points_ref , x_order, y_order);
410 interp_points_ref.resize(numInterpPoints,cellDim);
411
412 RealSpaceTools<double>::add(interpolant,exact_solution);
413
414 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
415
416 *outStream << "\nNorm-2 error between scalar components of exact solution polynomial of order ("
417 << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
418 << nrm << "\n";
419
420 if (nrm > zero) {
421 *outStream << "\n\nPatch test failed for solution polynomial order ("
422 << x_order << ", " << y_order << ") and basis order (scalar / vector) ("
423 << basis_order << ", " << basis_order + 1 << ")\n\n";
424 errorFlag++;
425 }
426
427 }
428 }
429 }
430
431 }
432 catch (const std::logic_error & err) {
433 *outStream << err.what() << "\n\n";
434 errorFlag = -1000;
435 };
436
437 if (errorFlag != 0)
438 std::cout << "End Result: TEST FAILED\n";
439 else
440 std::cout << "End Result: TEST PASSED\n";
441
442 // reset format state of std::cout
443 std::cout.copyfmt(oldFormatState);
444
445 return errorFlag;
446}
void u_exact(FieldContainer< double > &, const FieldContainer< double > &, int, int)
exact solution
Definition: test_02.cpp:91
void rhsFunc(FieldContainer< double > &, const FieldContainer< double > &, int, int)
right-hand side function
Definition: test_02.cpp:73
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::HDIV_QUAD_In_FEM class.
Header file for the Intrepid::HGRAD_QUAD_Cn_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 Quadrilateral cell.
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getReferenceSideNormal(ArraySideNormal &refSideNormal, const int &sideOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
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 Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.
static void add(Scalar *sumArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Adds contiguous data inArray1 and inArray2 of size size: sumArray = inArray1 + inArray2.