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_LINE_Cn_FEM.hpp"
53 #include "Intrepid_ArrayTools.hpp"
55 #include "Intrepid_CellTools.hpp"
56 #include "Intrepid_PointTools.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 
64 using namespace std;
65 using namespace Intrepid;
66 
67 void rhsFunc(FieldContainer<double> &, const FieldContainer<double> &, int);
68 void neumann(FieldContainer<double> &, const FieldContainer<double> &, const FieldContainer<double> &, int);
69 void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int);
70 
72 void 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 
98 void 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 
120 void 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 
131 int 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 }
Implementation of basic linear algebra functionality in Euclidean space.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Header file for the Intrepid::CellTools class.
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
int dimension(const int whichDim) const
Returns the specified dimension.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
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.
A stateless class for operations on cell data. Provides methods for: