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
50#include "Intrepid_HGRAD_LINE_Cn_FEM.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);
68void neumann(FieldContainer<double> &, const FieldContainer<double> &, const FieldContainer<double> &, int);
69void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int);
70
72void rhsFunc(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
73 if (degree == 0) {
74 for (int cell=0; cell<result.dimension(0); cell++) {
75 for (int pt=0; pt<result.dimension(1); pt++) {
76 result(cell,pt) = 1.0;
77
78 }
79 }
80 }
81 else if (degree == 1) {
82 for (int cell=0; cell<result.dimension(0); cell++) {
83 for (int pt=0; pt<result.dimension(1); pt++) {
84 result(cell,pt) = points(cell,pt,0);
85 }
86 }
87 }
88 else {
89 for (int cell=0; cell<result.dimension(0); cell++) {
90 for (int pt=0; pt<result.dimension(1); pt++) {
91 result(cell,pt) = pow(points(cell,pt,0), degree) - degree*(degree-1)*pow(points(cell,pt,0), degree-2);
92 }
93 }
94 }
95}
96
98void neumann(FieldContainer<double> & g_phi, const FieldContainer<double> & phi1, const FieldContainer<double> & phi2, int degree) {
99 double g_at_one, g_at_minus_one;
100 int num_fields;
101
102 if (degree == 0) {
103 g_at_one = 0.0;
104 g_at_minus_one = 0.0;
105 }
106 else {
107 g_at_one = degree*pow(1.0, degree-1);
108 g_at_minus_one = degree*pow(-1.0, degree-1);
109 }
110
111 num_fields = phi1.dimension(0);
112
113 for (int i=0; i<num_fields; i++) {
114 g_phi(0,i) = g_at_minus_one*phi1(i,0);
115 g_phi(1,i) = g_at_one*phi2(i,0);
116 }
117}
118
120void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
121 for (int cell=0; cell<result.dimension(0); cell++) {
122 for (int pt=0; pt<result.dimension(1); pt++) {
123 result(cell,pt) = pow(points(pt,0), degree);
124 }
125 }
126}
127
128
129
130
131int main(int argc, char *argv[]) {
132
133 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
134
135 // This little trick lets us print to std::cout only if
136 // a (dummy) command-line argument is provided.
137 int iprint = argc - 1;
138 Teuchos::RCP<std::ostream> outStream;
139 Teuchos::oblackholestream bhs; // outputs nothing
140 if (iprint > 0)
141 outStream = Teuchos::rcp(&std::cout, false);
142 else
143 outStream = Teuchos::rcp(&bhs, false);
144
145 // Save the format state of the original std::cout.
146 Teuchos::oblackholestream oldFormatState;
147 oldFormatState.copyfmt(std::cout);
148
149 *outStream \
150 << "===============================================================================\n" \
151 << "| |\n" \
152 << "| Unit Test (Basis_HGRAD_LINE_Cn_FEM) |\n" \
153 << "| |\n" \
154 << "| 1) Patch test involving mass and stiffness matrices, |\n" \
155 << "| for the Neumann problem on a REFERENCE line: |\n" \
156 << "| |\n" \
157 << "| - u'' + u = f in (-1,1), u' = g at -1,1 |\n" \
158 << "| |\n" \
159 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
160 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
161 << "| Robert Kirby (robert.c.kirby@ttu.gov), |\n" \
162 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
163 << "| |\n" \
164 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
165 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
166 << "| |\n" \
167 << "===============================================================================\n"\
168 << "| TEST 1: Patch test |\n"\
169 << "===============================================================================\n";
170
171
172 int errorFlag = 0;
173 double zero = 100*INTREPID_TOL;
174 outStream -> precision(20);
175
176
177 try {
178
179 int max_order = 10; // max total order of polynomial solution
180
181 // Define array containing points at which the solution is evaluated
182 int numIntervals = 100;
183 int numInterpPoints = numIntervals + 1;
184 FieldContainer<double> interp_points(numInterpPoints, 1);
185 for (int i=0; i<numInterpPoints; i++) {
186 interp_points(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
187 }
188
189 DefaultCubatureFactory<double> cubFactory; // create factory
190 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology
191
192
193 for (int basis_order=1; basis_order <= max_order; basis_order++) {
194 //create basis
195 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
196 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_SPECTRAL) );
197
198 for (int soln_order=1; soln_order <= basis_order; soln_order++) {
199 // evaluate exact solution
200 FieldContainer<double> exact_solution(1, numInterpPoints);
201 u_exact(exact_solution, interp_points, soln_order);
202
203 int numFields = lineBasis->getCardinality();
204
205 // create cubature
206 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, 2*basis_order-2);
207 int numCubPoints = lineCub->getNumPoints();
208 int spaceDim = lineCub->getDimension();
209
210 /* Computational arrays. */
211 FieldContainer<double> cub_points(numCubPoints, spaceDim);
212 FieldContainer<double> cub_points_physical(1, numCubPoints, spaceDim);
213 FieldContainer<double> cub_weights(numCubPoints);
214 FieldContainer<double> cell_nodes(1, 2, spaceDim);
215 FieldContainer<double> jacobian(1, numCubPoints, spaceDim, spaceDim);
216 FieldContainer<double> jacobian_inv(1, numCubPoints, spaceDim, spaceDim);
217 FieldContainer<double> jacobian_det(1, numCubPoints);
218 FieldContainer<double> weighted_measure(1, numCubPoints);
219
220 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints);
221 FieldContainer<double> transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
222 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
223 FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
224 FieldContainer<double> transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
225 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
226 FieldContainer<double> fe_matrix(1, numFields, numFields);
227
228 FieldContainer<double> rhs_at_cub_points_physical(1, numCubPoints);
229 FieldContainer<double> rhs_and_soln_vector(1, numFields);
230
231 FieldContainer<double> one_point(1, 1);
232 FieldContainer<double> value_of_basis_at_one(numFields, 1);
233 FieldContainer<double> value_of_basis_at_minusone(numFields, 1);
234 FieldContainer<double> bc_neumann(2, numFields);
235
236 FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
237 FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
238 FieldContainer<double> interpolant(1, numInterpPoints);
239
240 FieldContainer<int> ipiv(numFields);
241
242 /******************* START COMPUTATION ***********************/
243
244 // get cubature points and weights
245 lineCub->getCubature(cub_points, cub_weights);
246
247 // fill cell vertex array
248 cell_nodes(0, 0, 0) = -1.0;
249 cell_nodes(0, 1, 0) = 1.0;
250
251 // compute geometric cell information
252 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, line);
253 CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
254 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
255
256 // compute weighted measure
257 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
258
260 // Computing mass matrices:
261 // tabulate values of basis functions at (reference) cubature points
262 lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
263
264 // transform values of basis functions
265 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
266 value_of_basis_at_cub_points);
267
268 // multiply with weighted measure
269 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
270 weighted_measure,
271 transformed_value_of_basis_at_cub_points);
272
273 // compute mass matrices
274 FunctionSpaceTools::integrate<double>(fe_matrix,
275 transformed_value_of_basis_at_cub_points,
276 weighted_transformed_value_of_basis_at_cub_points,
277 COMP_CPP);
279
281 // Computing stiffness matrices:
282 // tabulate gradients of basis functions at (reference) cubature points
283 lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
284
285 // transform gradients of basis functions
286 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
287 jacobian_inv,
288 grad_of_basis_at_cub_points);
289
290 // multiply with weighted measure
291 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
292 weighted_measure,
293 transformed_grad_of_basis_at_cub_points);
294
295 // compute stiffness matrices and sum into fe_matrix
296 FunctionSpaceTools::integrate<double>(fe_matrix,
297 transformed_grad_of_basis_at_cub_points,
298 weighted_transformed_grad_of_basis_at_cub_points,
299 COMP_CPP,
300 true);
302
304 // Computing RHS contributions:
305 // map (reference) cubature points to physical space
306 CellTools<double>::mapToPhysicalFrame(cub_points_physical, cub_points, cell_nodes, line);
307
308 // evaluate rhs function
309 rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
310
311 // compute rhs
312 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
313 rhs_at_cub_points_physical,
314 weighted_transformed_value_of_basis_at_cub_points,
315 COMP_CPP);
316
317 // compute neumann b.c. contributions and adjust rhs
318 one_point(0,0) = 1.0; lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
319 one_point(0,0) = -1.0; lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
320 neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
321 for (int i=0; i<numFields; i++) {
322 rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
323 rhs_and_soln_vector(0, i) += bc_neumann(1, i);
324 }
326
328 // Solution of linear system:
329 int info = 0;
330 Teuchos::LAPACK<int, double> solver;
331 //solver.GESV(numRows, 1, &fe_mat(0,0), numRows, &ipiv(0), &fe_vec(0), numRows, &info);
332 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
334
336 // Building interpolant:
337 // evaluate basis at interpolation points
338 lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
339 // transform values of basis functions
340 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
341 value_of_basis_at_interp_points);
342 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
344
345 /******************* END COMPUTATION ***********************/
346
347 RealSpaceTools<double>::subtract(interpolant, exact_solution);
348
349 *outStream << "\nNorm-2 difference between exact solution polynomial of order "
350 << soln_order << " and finite element interpolant of order " << basis_order << ": "
351 << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) << "\n";
352
353 if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) > zero) {
354 *outStream << "\n\nPatch test failed for solution polynomial order "
355 << soln_order << " and basis order " << basis_order << "\n\n";
356 errorFlag++;
357 }
358
359 } // end for basis_order
360
361 } // end for soln_order
362
363 }
364 // Catch unexpected errors
365 catch (const std::logic_error & err) {
366 *outStream << err.what() << "\n\n";
367 errorFlag = -1000;
368 };
369
370 if (errorFlag != 0)
371 std::cout << "End Result: TEST FAILED\n";
372 else
373 std::cout << "End Result: TEST PASSED\n";
374
375 // reset format state of std::cout
376 std::cout.copyfmt(oldFormatState);
377
378 return errorFlag;
379}
void u_exact(FieldContainer< double > &, const FieldContainer< double > &, int)
exact solution
Definition test_02.cpp:120
void neumann(FieldContainer< double > &, const FieldContainer< double > &, const FieldContainer< double > &, int)
neumann boundary conditions
Definition test_02.cpp:98
void rhsFunc(FieldContainer< double > &, const FieldContainer< double > &, int)
right-hand side function
Definition test_02.cpp:72
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 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 locally H(grad)-compatible FEM basis of variable order on the [-1,...
A stateless class for operations on cell data. Provides methods for:
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
Implementation of basic linear algebra functionality in Euclidean space.
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.