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