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