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
50#include "Intrepid_HGRAD_TRI_Cn_FEM_ORTH.hpp"
57#include "Teuchos_oblackholestream.hpp"
58#include "Teuchos_RCP.hpp"
59#include "Teuchos_GlobalMPISession.hpp"
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_LAPACK.hpp"
63
64using namespace std;
65using namespace Intrepid;
66
67void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int );
68void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int, int );
69
70// This is the rhs for (div tau,w) = (f,w),
71// which makes f the negative Laplacian of scalar solution
72void rhsFunc( FieldContainer<double> &result,
73 const FieldContainer<double> &points,
74 int xd,
75 int yd )
76{
77 for (int cell=0;cell<result.dimension(0);cell++) {
78 for (int pt=0;pt<result.dimension(1);pt++) {
79 result(cell,pt) = 0.0;
80 if (xd >=2) {
81 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd);
82 }
83 if (yd >=2) {
84 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2);
85 }
86 }
87 }
88}
89
90void u_exact( FieldContainer<double> &result,
91 const FieldContainer<double> &points,
92 int xd,
93 int yd)
94{
95 for (int cell=0;cell<result.dimension(0);cell++){
96 for (int pt=0;pt<result.dimension(1);pt++) {
97 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd);
98 }
99 }
100 return;
101}
102
103int main(int argc, char *argv[]) {
104 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
105
106 // This little trick lets us print to std::cout only if
107 // a (dummy) command-line argument is provided.
108 int iprint = argc - 1;
109 Teuchos::RCP<std::ostream> outStream;
110 Teuchos::oblackholestream bhs; // outputs nothing
111 if (iprint > 0)
112 outStream = Teuchos::rcp(&std::cout, false);
113 else
114 outStream = Teuchos::rcp(&bhs, false);
115
116 // Save the format state of the original std::cout.
117 Teuchos::oblackholestream oldFormatState;
118 oldFormatState.copyfmt(std::cout);
119
120 *outStream \
121 << "===============================================================================\n" \
122 << "| |\n" \
123 << "| Unit Test (Basis_HGRAD_TRI_Cn_FEM_ORTH) |\n" \
124 << "| |\n" \
125 << "| 1) Patch test involving H(div) matrices |\n" \
126 << "| for the Dirichlet problem on a triangular patch |\n" \
127 << "| Omega with boundary Gamma. |\n" \
128 << "| |\n" \
129 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
130 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
131 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
132 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
133 << "| |\n" \
134 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
135 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
136 << "| |\n" \
137 << "===============================================================================\n"\
138 << "| TEST 1: Patch test |\n"\
139 << "===============================================================================\n";
140
141
142 int errorFlag = 0;
143
144 outStream -> precision(16);
145
146 try {
147 DefaultCubatureFactory<double> cubFactory; // create cubature factory
148 shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >()); // create parent cell topology
149 shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >()); // create relevant subcell (side) topology
150
151 int cellDim = cell.getDimension();
152 int sideDim = side.getDimension();
153
154 int min_order = 0;
155 int max_order = 9;
156
157 int numIntervals = max_order;
158 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2;
159 FieldContainer<double> interp_points_ref(numInterpPoints, 2);
160 int counter = 0;
161 for (int j=0; j<=numIntervals; j++) {
162 for (int i=0; i<=numIntervals; i++) {
163 if (i <= numIntervals-j) {
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
171 interp_points_ref.resize(numInterpPoints,2);
172
173 for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
174 // create bases
175 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
176 Teuchos::rcp(new Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >(basis_order+1,POINTTYPE_EQUISPACED) );
177 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
178 Teuchos::rcp(new Basis_HGRAD_TRI_Cn_FEM_ORTH<double,FieldContainer<double> >(basis_order) );
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-x_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<3;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}
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_TRI_In_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Implementation of the default H(div)-compatible Raviart-Thomas basis of arbitrary degree on Triangle ...
Implementation of the default H(grad)-compatible orthogonal basis (Dubiner) of arbitrary degree on tr...
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.
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.