Intrepid
example_08.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
81// Intrepid includes
89#include "Intrepid_Utils.hpp"
90
91// Epetra includes
92#include "Epetra_Time.h"
93#include "Epetra_Map.h"
94#include "Epetra_FECrsMatrix.h"
95#include "Epetra_FEVector.h"
96#include "Epetra_SerialComm.h"
97
98// Teuchos includes
99#include "Teuchos_oblackholestream.hpp"
100#include "Teuchos_RCP.hpp"
101#include "Teuchos_BLAS.hpp"
102
103// Shards includes
104#include "Shards_CellTopology.hpp"
105
106// EpetraExt includes
107#include "EpetraExt_RowMatrixOut.h"
108#include "EpetraExt_MultiVectorOut.h"
109
110using namespace std;
111using namespace Intrepid;
112
113// Functions to evaluate exact solution and derivatives
114double evalu(double & x, double & y, double & z);
115int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
116double evalDivGradu(double & x, double & y, double & z);
117
118int main(int argc, char *argv[]) {
119
120 //Check number of arguments
121 if (argc < 4) {
122 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
123 std::cout <<"Usage:\n\n";
124 std::cout <<" ./Intrepid_example_Drivers_Example_05.exe deg NX NY verbose\n\n";
125 std::cout <<" where \n";
126 std::cout <<" int deg - polynomial degree to be used (assumed > 1) \n";
127 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
128 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
129 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
130 exit(1);
131 }
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 > 2)
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 << "| Example: Generate Stiffness Matrix and Right Hand Side Vector for |\n" \
151 << "| Poisson Equation on Quadrilateral Mesh |\n" \
152 << "| |\n" \
153 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
154 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
155 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
156 << "| |\n" \
157 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
158 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
159 << "| |\n" \
160 << "===============================================================================\n";
161
162
163 // ************************************ GET INPUTS **************************************
164
165 int deg = atoi(argv[1]); // polynomial degree to use
166 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
167 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
168
169
170 // *********************************** CELL TOPOLOGY **********************************
171
172 // Get cell topology for base hexahedron
173 typedef shards::CellTopology CellTopology;
174 CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
175
176 // Get dimensions
177 int numNodesPerElem = quad_4.getNodeCount();
178 int spaceDim = quad_4.getDimension();
179
180 // *********************************** GENERATE MESH ************************************
181
182 *outStream << "Generating mesh ... \n\n";
183
184 *outStream << " NX" << " NY\n";
185 *outStream << std::setw(5) << NX <<
186 std::setw(5) << NY << "\n\n";
187
188 // Print mesh information
189 int numElems = NX*NY;
190 int numNodes = (NX+1)*(NY+1);
191 *outStream << " Number of Elements: " << numElems << " \n";
192 *outStream << " Number of Nodes: " << numNodes << " \n\n";
193
194 // Square
195 double leftX = 0.0, rightX = 1.0;
196 double leftY = 0.0, rightY = 1.0;
197
198 // Mesh spacing
199 double hx = (rightX-leftX)/((double)NX);
200 double hy = (rightY-leftY)/((double)NY);
201
202 // Get nodal coordinates
203 FieldContainer<double> nodeCoord(numNodes, spaceDim);
204 FieldContainer<int> nodeOnBoundary(numNodes);
205 int inode = 0;
206 for (int j=0; j<NY+1; j++) {
207 for (int i=0; i<NX+1; i++) {
208 nodeCoord(inode,0) = leftX + (double)i*hx;
209 nodeCoord(inode,1) = leftY + (double)j*hy;
210 if (j==0 || i==0 || j==NY || i==NX){
211 nodeOnBoundary(inode)=1;
212 }
213 else {
214 nodeOnBoundary(inode)=0;
215 }
216 inode++;
217 }
218 }
219#define DUMP_DATA
220#ifdef DUMP_DATA
221 // Print nodal coords
222 ofstream fcoordout("coords.dat");
223 for (int i=0; i<numNodes; i++) {
224 fcoordout << nodeCoord(i,0) <<" ";
225 fcoordout << nodeCoord(i,1) <<"\n";
226 }
227 fcoordout.close();
228#endif
229
230
231 // Element to Node map
232 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
233 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
234 int ielem=0;
235 for (int j=0; j<NY; j++) {
236 for (int i=0; i<NX; i++) {
237 elemToNode(ielem,0) = (NX + 1)*j + i;
238 elemToNode(ielem,1) = (NX + 1)*j + i + 1;
239 elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
240 elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
241 ielem++;
242 }
243 }
244#ifdef DUMP_DATA
245 // Output connectivity
246 ofstream fe2nout("elem2node.dat");
247 for (int j=0; j<NY; j++) {
248 for (int i=0; i<NX; i++) {
249 ielem = i + j * NX;
250 for (int m=0; m<numNodesPerElem; m++){
251 fe2nout << elemToNode(ielem,m) <<" ";
252 }
253 fe2nout <<"\n";
254 }
255 }
256 fe2nout.close();
257#endif
258
259
260 // ************************************ CUBATURE **************************************
261 *outStream << "Getting cubature ... \n\n";
262
263 // Get numerical integration points and weights
264 DefaultCubatureFactory<double> cubFactory;
265 int cubDegree = 2*deg;
266 Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree);
267
268 int cubDim = quadCub->getDimension();
269 int numCubPoints = quadCub->getNumPoints();
270
271 FieldContainer<double> cubPoints(numCubPoints, cubDim);
272 FieldContainer<double> cubWeights(numCubPoints);
273
274 quadCub->getCubature(cubPoints, cubWeights);
275
276
277 // ************************************** BASIS ***************************************
278
279 *outStream << "Getting basis ... \n\n";
280
281 // Define basis
282 Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
283 int numFieldsG = quadHGradBasis.getCardinality();
284 FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
285 FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
286
287 // Evaluate basis values and gradients at cubature points
288 quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
289 quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
290
291 // create the local-global mapping for higher order elements
292 FieldContainer<int> ltgMapping(numElems,numFieldsG);
293 const int numDOF = (NX*deg+1)*(NY*deg+1);
294 ielem=0;
295 for (int j=0;j<NY;j++) {
296 for (int i=0;i<NX;i++) {
297 const int start = deg * j * ( NX * deg + 1 ) + i * deg;
298 // loop over local dof on this cell
299 int local_dof_cur=0;
300 for (int vertical=0;vertical<=deg;vertical++) {
301 for (int horizontal=0;horizontal<=deg;horizontal++) {
302 ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
303 local_dof_cur++;
304 }
305 }
306 ielem++;
307 }
308 }
309#ifdef DUMP_DATA
310 // Output ltg mapping
311 ofstream ltgout("ltg.dat");
312 for (int j=0; j<NY; j++) {
313 for (int i=0; i<NX; i++) {
314 ielem = i + j * NX;
315 for (int m=0; m<numFieldsG; m++){
316 ltgout << ltgMapping(ielem,m) <<" ";
317 }
318 ltgout <<"\n";
319 }
320 }
321 ltgout.close();
322#endif
323
324 // ******** CREATE ALL LOCAL STIFFNESS MATRICES *********
325 *outStream << "Building stiffness matrix and right hand side ... \n\n";
326
327 // Settings and data structures for mass and stiffness matrices
328 typedef CellTools<double> CellTools;
329 typedef FunctionSpaceTools fst;
330 int numCells = numElems;
331
332 // Container for nodes
333 FieldContainer<double> cellVertices(numCells,numNodesPerElem,spaceDim);
334
335 // Containers for Jacobian
336 FieldContainer<double> cellJacobian(numCells, numCubPoints, spaceDim, spaceDim);
337 FieldContainer<double> cellJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
338 FieldContainer<double> cellJacobDet(numCells, numCubPoints);
339
340 // Containers for element HGRAD stiffness matrices
341 FieldContainer<double> localStiffMatrices(numCells, numFieldsG, numFieldsG);
342 FieldContainer<double> transformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
343 FieldContainer<double> weightedTransformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
344 FieldContainer<double> weightedMeasure(numCells, numCubPoints);
345
346
347 // Global arrays in Epetra format
348 Epetra_SerialComm Comm;
349 Epetra_Map globalMapG(numDOF, 0, Comm);
350
351 Epetra_Time graphTimer(Comm);
352 Epetra_CrsGraph grph( Copy , globalMapG , 4 * numFieldsG );
353 for (int k=0;k<numElems;k++)
354 {
355 for (int i=0;i<numFieldsG;i++)
356 {
357 grph.InsertGlobalIndices(ltgMapping(k,i),numFieldsG,&ltgMapping(k,0));
358 }
359 }
360 grph.FillComplete();
361
362 const double graphTime = graphTimer.ElapsedTime();
363 std::cout << "Graph computed in " << graphTime << "\n";
364
365 Epetra_Time instantiateTimer( Comm );
366 Epetra_FECrsMatrix StiffMatrix( Copy , grph );
367 const double instantiateTime = instantiateTimer.ElapsedTime( );
368 std::cout << "Matrix instantiated in " << instantiateTime << "\n";
369
370 Epetra_FEVector u(globalMapG);
371 Epetra_FEVector Ku(globalMapG);
372
373 u.Random();
374
375 // ************************** Compute element HGrad stiffness matrices *******************************
376 // Get vertices of all the cells
377
378 for (int i=0;i<numElems;i++)
379 {
380 for (int j=0;j<4;j++)
381 {
382 const int nodeCur = elemToNode(i,j);
383 for (int k=0;k<spaceDim;k++)
384 {
385 cellVertices(i,j,k) = nodeCoord(nodeCur,k);
386 }
387 }
388 }
389
390 Epetra_Time localConstructTimer(Comm);
391
392 // Compute all cell Jacobians, their inverses and their determinants
393 CellTools::setJacobian(cellJacobian, cubPoints, cellVertices, quad_4);
394 CellTools::setJacobianInv(cellJacobInv, cellJacobian );
395 CellTools::setJacobianDet(cellJacobDet, cellJacobian );
396
397 // transform reference element gradients to each cell
398 fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads);
399
400 // compute weighted measure
401 fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights);
402
403 // multiply values with weighted measure
404 fst::multiplyMeasure<double>(weightedTransformedBasisGradients,
405 weightedMeasure, transformedBasisGradients);
406
407 // integrate to compute element stiffness matrix
408 fst::integrate<double>(localStiffMatrices,
409 transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS);
410
411 const double localConstructTime = localConstructTimer.ElapsedTime();
412 std::cout << "Time to build local matrices (including Jacobian computation): "<< localConstructTime << "\n";
413
414 Epetra_Time insertionTimer(Comm);
415
416 // *** Element loop ***
417 for (int k=0; k<numElems; k++)
418 {
419 // assemble into global matrix
420 StiffMatrix.InsertGlobalValues(numFieldsG,&ltgMapping(k,0),numFieldsG,&ltgMapping(k,0),&localStiffMatrices(k,0,0));
421
422 }
423 StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
424 const double insertionTime = insertionTimer.ElapsedTime( );
425
426 std::cout << "Time to assemble global matrix from local matrices: " << insertionTime << "\n";
427
428
429#ifdef DUMP_DATA
430 // Dump matrices to disk
431// EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
432// EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
433#endif
434
435
436 std::cout << "End Result: TEST PASSED\n";
437
438 // reset format state of std::cout
439 std::cout.copyfmt(oldFormatState);
440
441 return 0;
442}
443
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.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...