Intrepid
example_05.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
83// Intrepid includes
91#include "Intrepid_Utils.hpp"
92
93// Epetra includes
94#include "Epetra_Time.h"
95#include "Epetra_Map.h"
96#include "Epetra_FECrsMatrix.h"
97#include "Epetra_FEVector.h"
98#include "Epetra_SerialComm.h"
99
100// Teuchos includes
101#include "Teuchos_oblackholestream.hpp"
102#include "Teuchos_RCP.hpp"
103#include "Teuchos_BLAS.hpp"
104
105// Shards includes
106#include "Shards_CellTopology.hpp"
107
108// EpetraExt includes
109#include "EpetraExt_RowMatrixOut.h"
110#include "EpetraExt_MultiVectorOut.h"
111
112using namespace std;
113using namespace Intrepid;
114
115// Functions to evaluate exact solution and derivatives
116double evalu(double & x, double & y, double & z);
117int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
118double evalDivGradu(double & x, double & y, double & z);
119
120int main(int argc, char *argv[]) {
121
122 //Check number of arguments
123 if (argc < 4) {
124 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
125 std::cout <<"Usage:\n\n";
126 std::cout <<" ./Intrepid_example_Drivers_Example_05.exe deg NX NY verbose\n\n";
127 std::cout <<" where \n";
128 std::cout <<" int deg - polynomial degree to be used (assumed > 1) \n";
129 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
130 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
131 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
132 exit(1);
133 }
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 > 2)
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 << "| Example: Generate Stiffness Matrix and Right Hand Side Vector for |\n" \
153 << "| Poisson Equation on Quadrilateral Mesh |\n" \
154 << "| |\n" \
155 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
156 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
157 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
158 << "| |\n" \
159 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
160 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
161 << "| |\n" \
162 << "===============================================================================\n";
163
164
165 // ************************************ GET INPUTS **************************************
166
167 int deg = atoi(argv[1]); // polynomial degree to use
168 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
169 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
170
171
172 // *********************************** CELL TOPOLOGY **********************************
173
174 // Get cell topology for base hexahedron
175 typedef shards::CellTopology CellTopology;
176 CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
177
178 // Get dimensions
179 int numNodesPerElem = quad_4.getNodeCount();
180 int spaceDim = quad_4.getDimension();
181
182 // *********************************** GENERATE MESH ************************************
183
184 *outStream << "Generating mesh ... \n\n";
185
186 *outStream << " NX" << " NY\n";
187 *outStream << std::setw(5) << NX <<
188 std::setw(5) << NY << "\n\n";
189
190 // Print mesh information
191 int numElems = NX*NY;
192 int numNodes = (NX+1)*(NY+1);
193 *outStream << " Number of Elements: " << numElems << " \n";
194 *outStream << " Number of Nodes: " << numNodes << " \n\n";
195
196 // Square
197 double leftX = 0.0, rightX = 1.0;
198 double leftY = 0.0, rightY = 1.0;
199
200 // Mesh spacing
201 double hx = (rightX-leftX)/((double)NX);
202 double hy = (rightY-leftY)/((double)NY);
203
204 // Get nodal coordinates
205 FieldContainer<double> nodeCoord(numNodes, spaceDim);
206 FieldContainer<int> nodeOnBoundary(numNodes);
207 int inode = 0;
208 for (int j=0; j<NY+1; j++) {
209 for (int i=0; i<NX+1; i++) {
210 nodeCoord(inode,0) = leftX + (double)i*hx;
211 nodeCoord(inode,1) = leftY + (double)j*hy;
212 if (j==0 || i==0 || j==NY || i==NX){
213 nodeOnBoundary(inode)=1;
214 }
215 else {
216 nodeOnBoundary(inode)=0;
217 }
218 inode++;
219 }
220 }
221#define DUMP_DATA
222#ifdef DUMP_DATA
223 // Print nodal coords
224 ofstream fcoordout("coords.dat");
225 for (int i=0; i<numNodes; i++) {
226 fcoordout << nodeCoord(i,0) <<" ";
227 fcoordout << nodeCoord(i,1) <<"\n";
228 }
229 fcoordout.close();
230#endif
231
232
233 // Element to Node map
234 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
235 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
236 int ielem = 0;
237 for (int j=0; j<NY; j++) {
238 for (int i=0; i<NX; i++) {
239 elemToNode(ielem,0) = (NX + 1)*j + i;
240 elemToNode(ielem,1) = (NX + 1)*j + i + 1;
241 elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
242 elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
243 ielem++;
244 }
245 }
246#ifdef DUMP_DATA
247 // Output connectivity
248 ofstream fe2nout("elem2node.dat");
249 for (int j=0; j<NY; j++) {
250 for (int i=0; i<NX; i++) {
251 ielem = i + j * NX;
252 for (int m=0; m<numNodesPerElem; m++){
253 fe2nout << elemToNode(ielem,m) <<" ";
254 }
255 fe2nout <<"\n";
256 }
257 }
258 fe2nout.close();
259#endif
260
261
262 // ************************************ CUBATURE **************************************
263 *outStream << "Getting cubature ... \n\n";
264
265 // Get numerical integration points and weights
267 int cubDegree = 2*deg;
268 Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree);
269
270 int cubDim = quadCub->getDimension();
271 int numCubPoints = quadCub->getNumPoints();
272
273 FieldContainer<double> cubPoints(numCubPoints, cubDim);
274 FieldContainer<double> cubWeights(numCubPoints);
275
276 quadCub->getCubature(cubPoints, cubWeights);
277
278
279 // ************************************** BASIS ***************************************
280
281 *outStream << "Getting basis ... \n\n";
282
283 // Define basis
284 Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
285 int numFieldsG = quadHGradBasis.getCardinality();
286 FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
287 FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
288
289 // Evaluate basis values and gradients at cubature points
290 quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
291 quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
292
293 // create the local-global mapping for higher order elements
294 FieldContainer<int> ltgMapping(numElems,numFieldsG);
295 const int numDOF = (NX*deg+1)*(NY*deg+1);
296 ielem=0;
297 for (int j=0;j<NY;j++) {
298 for (int i=0;i<NX;i++) {
299 const int start = deg * j * ( NX * deg + 1 ) + i * deg;
300 // loop over local dof on this cell
301 int local_dof_cur=0;
302 for (int vertical=0;vertical<=deg;vertical++) {
303 for (int horizontal=0;horizontal<=deg;horizontal++) {
304 ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
305 local_dof_cur++;
306 }
307 }
308 ielem++;
309 }
310 }
311#ifdef DUMP_DATA
312 // Output ltg mapping
313 ofstream ltgout("ltg.dat");
314 for (int j=0; j<NY; j++) {
315 for (int i=0; i<NX; i++) {
316 ielem = i + j * NX;
317 for (int m=0; m<numFieldsG; m++){
318 ltgout << ltgMapping(ielem,m) <<" ";
319 }
320 ltgout <<"\n";
321 }
322 }
323 ltgout.close();
324#endif
325
326 // ******** CREATE A SINGLE STIFFNESS MATRIX, WHICH IS REPLICATED ON ALL ELEMENTS *********
327 *outStream << "Building stiffness matrix and right hand side ... \n\n";
328
329 // Settings and data structures for mass and stiffness matrices
331 typedef FunctionSpaceTools fst;
332 int numCells = 1;
333
334 // Container for nodes
335 FieldContainer<double> refQuadNodes(numCells, numNodesPerElem, spaceDim);
336 // Containers for Jacobian
337 FieldContainer<double> refQuadJacobian(numCells, numCubPoints, spaceDim, spaceDim);
338 FieldContainer<double> refQuadJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
339 FieldContainer<double> refQuadJacobDet(numCells, numCubPoints);
340 // Containers for element HGRAD stiffness matrix
341 FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
342 FieldContainer<double> weightedMeasure(numCells, numCubPoints);
343 FieldContainer<double> quadGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
344 FieldContainer<double> quadGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
345 // Containers for right hand side vectors
346 FieldContainer<double> rhsData(numCells, numCubPoints);
347 FieldContainer<double> localRHS(numCells, numFieldsG);
348 FieldContainer<double> quadGValsTransformed(numCells, numFieldsG, numCubPoints);
349 FieldContainer<double> quadGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
350 // Container for cubature points in physical space
351 FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
352
353 // Global arrays in Epetra format
354 Epetra_SerialComm Comm;
355 Epetra_Map globalMapG(numDOF, 0, Comm);
356 Epetra_Time instantiateTimer(Comm);
357 Epetra_FECrsMatrix StiffMatrix(Copy, globalMapG, 4*numFieldsG);
358 const double instantiateTime = instantiateTimer.ElapsedTime();
359 std::cout << "Time to instantiate sparse matrix " << instantiateTime << "\n";
360 Epetra_FEVector u(globalMapG);
361 Epetra_FEVector Ku(globalMapG);
362
363 u.Random();
364
365 // ************************** Compute element HGrad stiffness matrices *******************************
366 refQuadNodes(0,0,0) = 0.0;
367 refQuadNodes(0,0,1) = 0.0;
368 refQuadNodes(0,1,0) = hx;
369 refQuadNodes(0,1,1) = 0.0;
370 refQuadNodes(0,2,0) = hx;
371 refQuadNodes(0,2,1) = hy;
372 refQuadNodes(0,3,0) = 0.0;
373 refQuadNodes(0,3,1) = hy;
374
375 // Compute cell Jacobians, their inverses and their determinants
376 CellTools::setJacobian(refQuadJacobian, cubPoints, refQuadNodes, quad_4);
377 CellTools::setJacobianInv(refQuadJacobInv, refQuadJacobian );
378 CellTools::setJacobianDet(refQuadJacobDet, refQuadJacobian );
379
380 // transform from [-1,1]^2 to [0,hx]x[0,hy]
381 fst::HGRADtransformGRAD<double>(quadGradsTransformed, refQuadJacobInv, quadGrads);
382
383 // compute weighted measure
384 fst::computeCellMeasure<double>(weightedMeasure, refQuadJacobDet, cubWeights);
385
386 // multiply values with weighted measure
387 fst::multiplyMeasure<double>(quadGradsTransformedWeighted,
388 weightedMeasure, quadGradsTransformed);
389
390 // integrate to compute element stiffness matrix
391 fst::integrate<double>(localStiffMatrix,
392 quadGradsTransformed, quadGradsTransformedWeighted, COMP_BLAS);
393
394
395 Epetra_Time assemblyTimer(Comm);
396
397 // *** Element loop ***
398 for (int k=0; k<numElems; k++)
399 {
400 // assemble into global matrix
401 StiffMatrix.InsertGlobalValues(numFieldsG,&ltgMapping(k,0),numFieldsG,&ltgMapping(k,0),&localStiffMatrix(0,0,0));
402
403 }
404
405
406 // Assemble global matrices
407 StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
408
409 double assembleTime = assemblyTimer.ElapsedTime();
410 std::cout << "Time to insert reference element matrix into global matrix: " << assembleTime << std::endl;
411 std::cout << "There are " << StiffMatrix.NumGlobalNonzeros() << " nonzeros in the matrix.\n";
412 std::cout << "There are " << numDOF << " global degrees of freedom.\n";
413
414 Epetra_Time multTimer(Comm);
415 StiffMatrix.Apply(u,Ku);
416 double multTime = multTimer.ElapsedTime();
417 std::cout << "Time to apply: " << multTime << std::endl;
418
419// // Adjust stiffness matrix and rhs based on boundary conditions
420// for (int row = 0; row<numNodes; row++){
421// if (nodeOnBoundary(row)) {
422// int rowindex = row;
423// for (int col=0; col<numNodes; col++){
424// double val = 0.0;
425// int colindex = col;
426// StiffMatrix.ReplaceGlobalValues(1, &rowindex, 1, &colindex, &val);
427// }
428// double val = 1.0;
429// StiffMatrix.ReplaceGlobalValues(1, &rowindex, 1, &rowindex, &val);
430// val = 0.0;
431// rhs.ReplaceGlobalValues(1, &rowindex, &val);
432// }
433// }
434
435#ifdef DUMP_DATA
436 // Dump matrices to disk
437// EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
438// EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
439#endif
440
441 std::cout << "End Result: TEST PASSED\n";
442
443 // reset format state of std::cout
444 std::cout.copyfmt(oldFormatState);
445
446 return 0;
447}
448
449
450// Calculates value of exact solution u
451 double evalu(double & x, double & y, double & z)
452 {
453 /*
454 // function1
455 double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
456 */
457
458 // function2
459 double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
460
461 return exactu;
462 }
463
464// Calculates gradient of exact solution u
465 int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3)
466 {
467 /*
468 // function 1
469 gradu1 = M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
470 gradu2 = M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
471 gradu3 = M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
472 */
473
474 // function2
475 gradu1 = (M_PI*cos(M_PI*x)+sin(M_PI*x))
476 *sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
477 gradu2 = (M_PI*cos(M_PI*y)+sin(M_PI*y))
478 *sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z);
479 gradu3 = (M_PI*cos(M_PI*z)+sin(M_PI*z))
480 *sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z);
481
482 return 0;
483 }
484
485// Calculates Laplacian of exact solution u
486 double evalDivGradu(double & x, double & y, double & z)
487 {
488 /*
489 // function 1
490 double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
491 */
492
493 // function 2
494 double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
495 + 2.0*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
496 + 2.0*M_PI*cos(M_PI*y)*sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z)
497 + 2.0*M_PI*cos(M_PI*z)*sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z)
498 + 3.0*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
499
500 return divGradu;
501 }
502
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::HGRAD_QUAD_Cn_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid utilities.
A stateless class for operations on cell data. Provides methods for:
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.
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....
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...