Intrepid
example_09.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
81// 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
114int main(int argc, char *argv[]) {
115
116 //Check number of arguments
117 if (argc < 4) {
118 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
119 std::cout <<"Usage:\n\n";
120 std::cout <<" ./Intrepid_example_Drivers_Example_09.exe deg NX NY verbose\n\n";
121 std::cout <<" where \n";
122 std::cout <<" int deg - polynomial degree to be used (assumed > 1) \n";
123 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
124 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
125 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
126 exit(1);
127 }
128
129 // This little trick lets us print to std::cout only if
130 // a (dummy) command-line argument is provided.
131 int iprint = argc - 1;
132 Teuchos::RCP<std::ostream> outStream;
133 Teuchos::oblackholestream bhs; // outputs nothing
134 if (iprint > 2)
135 outStream = Teuchos::rcp(&std::cout, false);
136 else
137 outStream = Teuchos::rcp(&bhs, false);
138
139 // Save the format state of the original std::cout.
140 Teuchos::oblackholestream oldFormatState;
141 oldFormatState.copyfmt(std::cout);
142
143 *outStream \
144 << "===============================================================================\n" \
145 << "| |\n" \
146 << "| Example: Apply Stiffness Matrix for |\n" \
147 << "| Poisson Equation on Quadrilateral Mesh |\n" \
148 << "| |\n" \
149 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
150 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
151 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
152 << "| |\n" \
153 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
154 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
155 << "| |\n" \
156 << "===============================================================================\n";
157
158
159 // ************************************ GET INPUTS **************************************
160
161 int deg = atoi(argv[1]); // polynomial degree to use
162 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
163 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
164
165
166 // *********************************** CELL TOPOLOGY **********************************
167
168 // Get cell topology for base hexahedron
169 typedef shards::CellTopology CellTopology;
170 CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
171
172 // Get dimensions
173 int numNodesPerElem = quad_4.getNodeCount();
174 int spaceDim = quad_4.getDimension();
175
176 // *********************************** GENERATE MESH ************************************
177
178 *outStream << "Generating mesh ... \n\n";
179
180 *outStream << " NX" << " NY\n";
181 *outStream << std::setw(5) << NX <<
182 std::setw(5) << NY << "\n\n";
183
184 // Print mesh information
185 int numElems = NX*NY;
186 int numNodes = (NX+1)*(NY+1);
187 *outStream << " Number of Elements: " << numElems << " \n";
188 *outStream << " Number of Nodes: " << numNodes << " \n\n";
189
190 // Square
191 double leftX = 0.0, rightX = 1.0;
192 double leftY = 0.0, rightY = 1.0;
193
194 // Mesh spacing
195 double hx = (rightX-leftX)/((double)NX);
196 double hy = (rightY-leftY)/((double)NY);
197
198 // Get nodal coordinates
199 FieldContainer<double> nodeCoord(numNodes, spaceDim);
200 FieldContainer<int> nodeOnBoundary(numNodes);
201 int inode = 0;
202 for (int j=0; j<NY+1; j++) {
203 for (int i=0; i<NX+1; i++) {
204 nodeCoord(inode,0) = leftX + (double)i*hx;
205 nodeCoord(inode,1) = leftY + (double)j*hy;
206 if (j==0 || i==0 || j==NY || i==NX){
207 nodeOnBoundary(inode)=1;
208 }
209 else {
210 nodeOnBoundary(inode)=0;
211 }
212 inode++;
213 }
214 }
215#define DUMP_DATA
216#ifdef DUMP_DATA
217 // Print nodal coords
218 ofstream fcoordout("coords.dat");
219 for (int i=0; i<numNodes; i++) {
220 fcoordout << nodeCoord(i,0) <<" ";
221 fcoordout << nodeCoord(i,1) <<"\n";
222 }
223 fcoordout.close();
224#endif
225
226
227 // Element to Node map
228 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
229 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
230 int ielem = 0;
231 for (int j=0; j<NY; j++) {
232 for (int i=0; i<NX; i++) {
233 elemToNode(ielem,0) = (NX + 1)*j + i;
234 elemToNode(ielem,1) = (NX + 1)*j + i + 1;
235 elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
236 elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
237 ielem++;
238 }
239 }
240#ifdef DUMP_DATA
241 // Output connectivity
242 ofstream fe2nout("elem2node.dat");
243 for (int j=0; j<NY; j++) {
244 for (int i=0; i<NX; i++) {
245 ielem = i + j * NX;
246 for (int m=0; m<numNodesPerElem; m++){
247 fe2nout << elemToNode(ielem,m) <<" ";
248 }
249 fe2nout <<"\n";
250 }
251 }
252 fe2nout.close();
253#endif
254
255
256 // ************************************ CUBATURE **************************************
257 *outStream << "Getting cubature ... \n\n";
258
259 // Get numerical integration points and weights
260 // I only need this on the line since I'm doing tensor products
261 DefaultCubatureFactory<double> cubFactory;
262
263 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
264 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
265
266 const int numCubPoints = glcub->getNumPoints();
267
268 FieldContainer<double> cubPoints1D(numCubPoints, 1);
269 FieldContainer<double> cubWeights1D(numCubPoints);
270
271 glcub->getCubature(cubPoints1D,cubWeights1D);
272
273
274 // ************************************** BASIS ***************************************
275 *outStream << "Getting basis ... \n\n";
276
277 // Define basis: I only need this on the line also
278 Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineHGradBasis(deg,POINTTYPE_SPECTRAL);
279 int numLineFieldsG = lineHGradBasis.getCardinality();
280 FieldContainer<double> lineGrads(numLineFieldsG, numCubPoints, 1);
281
282 // Evaluate basis values and gradients at cubature points
283 lineHGradBasis.getValues(lineGrads, cubPoints1D, OPERATOR_GRAD);
284
285 // ************************************** LTG mapping **********************************
286 FieldContainer<int> ltgMapping(numElems,numLineFieldsG*numLineFieldsG);
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<numLineFieldsG; m++){
310 ltgout << ltgMapping(ielem,m) <<" ";
311 }
312 ltgout <<"\n";
313 }
314 }
315 ltgout.close();
316#endif
317
318
319 // Global arrays in Epetra format
320 Epetra_SerialComm Comm;
321 Epetra_Map globalMapG(numDOF, 0, Comm);
322
323 Epetra_FEVector u(globalMapG);
324 Epetra_FEVector Ku(globalMapG);
325
326 u.Random();
327
328
329 // ************************** Compute element HGrad stiffness matrices *******************************
330// // Get vertices of all the cells
331// for (int i=0;i<numElems;i++)
332// {
333// for (int j=0;j<4;j++)
334// {
335// const int nodeCur = elemToNode(i,j);
336// for (int k=0;k<spaceDim;k++)
337// {
338// cellVertices(i,j,k) = nodeCoord(nodeCur,k);
339// }
340// }
341// }
342
343 FieldContainer<double> uScattered(numElems,numLineFieldsG*numLineFieldsG);
344 FieldContainer<double> KuScattered(numElems,numLineFieldsG*numLineFieldsG);
345
346 // need storage for derivatives of u on each cell
347 // the number of line dof should be the same as the
348 // number of cub points.
349 // This is indexed by Du(q2,q1)
350 FieldContainer<double> Du(numCubPoints,numCubPoints);
351
352
353
354 double *uVals = u[0];
355 double *KuVals = Ku[0];
356 Epetra_Time scatterTime(Comm);
357 *outStream << "Scattering\n";
358 // Scatter
359 for (int k=0; k<numElems; k++)
360 {
361 for (int i=0;i<numLineFieldsG*numLineFieldsG;i++)
362 {
363 uScattered(k,i) = uVals[ltgMapping(k,i)];
364 }
365 }
366 const double scatTime = scatterTime.ElapsedTime();
367 *outStream << "Scattered in time " << scatTime << "\n";
368
369 Epetra_Time applyTimer(Comm);
370
371 uScattered.resize(numElems,numLineFieldsG,numLineFieldsG);
372
373 for (int k=0;k<numElems;k++)
374 {
375 // local operation: x-derivative term of stiffness matrix
376 // evaluate x derivative of u on cell k.
377 for (int j=0;j<numLineFieldsG;j++)
378 {
379 for (int i=0;i<numLineFieldsG;i++)
380 {
381 Du(j,i) = 0.0;
382 for (int q=0;q<numLineFieldsG;q++)
383 {
384 Du(j,i) += uScattered(k,j,i) * lineGrads(i,q,0);
385 }
386 }
387 }
388
389 // initialize Ku
390 for (int i=0;i<numLineFieldsG*numLineFieldsG;i++)
391 {
392 KuScattered(k,i) = 0.0;
393 }
394
395 // loop over basis functions for x term
396 int cur = 0;
397 for (int j=0;j<numLineFieldsG;j++)
398 {
399 for (int i=0;i<numLineFieldsG;i++)
400 {
401 // do the quadrature
402 for (int q1=0;q1<numCubPoints;q1++)
403 {
404 KuScattered(k,cur) += cubWeights1D(j) * cubWeights1D(q1) * Du(j,q1) * lineGrads(i,q1,0);
405 }
406 cur ++;
407 }
408 }
409
410 // local operation: y-derivative term of stiffness matrix, store in Du
411 for (int j=0;j<numLineFieldsG;j++)
412 {
413 for (int i=0;i<numLineFieldsG;i++)
414 {
415 Du(j,i) = 0.0;
416 for (int q=0;q<numLineFieldsG;q++)
417 {
418 Du(j,i) += uScattered(k,j,i) * lineGrads(j,q,0);
419 }
420 }
421 }
422
423
424 // evaluate y-derivatives of u
425 cur = 0;
426 for (int j=0;j<numLineFieldsG;j++)
427 {
428 for (int i=0;i<numLineFieldsG;i++)
429 {
430 for (int q2=0;q2<numCubPoints;q2++)
431 {
432 KuScattered(k,cur) += cubWeights1D(q2) * Du(q2,i) * lineGrads(j,q2,0) * cubWeights1D(i);
433 }
434 }
435 }
436 }
437
438 uScattered.resize(numElems,numLineFieldsG*numLineFieldsG);
439
440 const double applyTime = applyTimer.ElapsedTime();
441
442 *outStream << "Local application: " << applyTime << "\n";
443
444 // gather
445 Epetra_Time gatherTimer(Comm);
446 // Gather
447 for (int k=0;k<numElems;k++)
448 {
449 for (int i=0;i<numLineFieldsG*numLineFieldsG;i++)
450 {
451 KuVals[ltgMapping(k,i)] += KuScattered(k,i);
452 }
453 }
454
455 const double gatherTime = gatherTimer.ElapsedTime();
456 *outStream << "Gathered in " << gatherTime << "\n";
457
458
459
460#ifdef DUMP_DATA
461 // Dump matrices to disk
462// EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
463// EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
464#endif
465
466
467 std::cout << "End Result: TEST PASSED\n";
468
469 // reset format state of std::cout
470 std::cout.copyfmt(oldFormatState);
471
472 return 0;
473}
474
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the Intrepid::CubaturePolylib 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.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...