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