88#include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
93#include "Epetra_Time.h"
94#include "Epetra_Map.h"
95#include "Epetra_FEVector.h"
96#include "Epetra_SerialComm.h"
99#include "Teuchos_oblackholestream.hpp"
100#include "Teuchos_RCP.hpp"
105#include "Shards_CellTopology.hpp"
108#include "EpetraExt_MultiVectorOut.h"
111using namespace Intrepid;
113int main(
int argc,
char *argv[]) {
117 std::cout <<
"\n>>> ERROR: Invalid number of arguments.\n\n";
118 std::cout <<
"Usage:\n\n";
119 std::cout <<
" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ 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 <<
" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
125 std::cout <<
" verbose (optional) - any character, indicates verbose output \n\n";
131 int iprint = argc - 1;
132 Teuchos::RCP<std::ostream> outStream;
133 Teuchos::oblackholestream bhs;
135 outStream = Teuchos::rcp(&std::cout,
false);
137 outStream = Teuchos::rcp(&bhs,
false);
140 Teuchos::oblackholestream oldFormatState;
141 oldFormatState.copyfmt(std::cout);
144 <<
"===============================================================================\n" \
146 <<
"| Example: Apply Stiffness Matrix for |\n" \
147 <<
"| Poisson Equation on Hexahedral Mesh |\n" \
149 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
150 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
151 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
153 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
154 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
156 <<
"===============================================================================\n";
161 int deg = atoi(argv[1]);
162 int NX = atoi(argv[2]);
163 int NY = atoi(argv[3]);
164 int NZ = atoi(argv[4]);
170 typedef shards::CellTopology CellTopology;
171 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
174 int numNodesPerElem = hex_8.getNodeCount();
175 int spaceDim = hex_8.getDimension();
179 *outStream <<
"Generating mesh ... \n\n";
181 *outStream <<
" NX" <<
" NY" <<
" NZ\n";
182 *outStream << std::setw(5) << NX <<
183 std::setw(5) << NY << std::setw(5) << NZ <<
"\n\n";
186 int numElems = NX*NY*NZ;
187 int numNodes = (NX+1)*(NY+1)*(NZ+1);
188 *outStream <<
" Number of Elements: " << numElems <<
" \n";
189 *outStream <<
" Number of Nodes: " << numNodes <<
" \n\n";
192 double leftX = 0.0, rightX = 1.0;
193 double leftY = 0.0, rightY = 1.0;
194 double leftZ = 0.0, rightZ = 1.0;
197 double hx = (rightX-leftX)/((
double)NX);
198 double hy = (rightY-leftY)/((
double)NY);
199 double hz = (rightZ-leftZ)/((
double)NZ);
205 for (
int k=0; k<NZ+1; k++)
207 for (
int j=0; j<NY+1; j++)
209 for (
int i=0; i<NX+1; i++)
211 nodeCoord(inode,0) = leftX + (double)i*hx;
212 nodeCoord(inode,1) = leftY + (double)j*hy;
213 nodeCoord(inode,2) = leftZ + (double)k*hz;
214 if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
216 nodeOnBoundary(inode)=1;
220 nodeOnBoundary(inode)=0;
229 ofstream fcoordout(
"coords.dat");
230 for (
int i=0; i<numNodes; i++) {
231 fcoordout << nodeCoord(i,0) <<
" ";
232 fcoordout << nodeCoord(i,1) <<
" ";
233 fcoordout << nodeCoord(i,2) <<
"\n";
243 for (
int k=0; k<NZ; k++)
245 for (
int j=0; j<NY; j++)
247 for (
int i=0; i<NX; i++)
249 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
250 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
251 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
252 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
253 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
254 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
255 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
256 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
263 ofstream fe2nout(
"elem2node.dat");
264 for (
int k=0;k<NZ;k++)
266 for (
int j=0; j<NY; j++)
268 for (
int i=0; i<NX; i++)
270 ielem = i + j * NX + k * NY * NY;
271 for (
int m=0; m<numNodesPerElem; m++)
273 fe2nout << elemToNode(ielem,m) <<
" ";
283 *outStream <<
"Getting cubature and basis ... \n\n";
290 const int numCubPoints = glcub->getNumPoints();
295 glcub->getCubature(cubPoints1D,cubWeights1D);
298 int numLineFieldsG = lineHGradBasis.getCardinality();
302 lineHGradBasis.getValues(lineGrads, cubPoints1D, OPERATOR_GRAD);
307 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
309 for (
int k=0;k<NZ;k++)
311 for (
int j=0;j<NY;j++)
313 for (
int i=0;i<NX;i++)
315 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
318 for (
int kloc=0;kloc<=deg;kloc++)
320 for (
int jloc=0;jloc<=deg;jloc++)
322 for (
int iloc=0;iloc<=deg;iloc++)
324 ltgMapping(ielem,local_dof_cur) = start
325 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
326 + jloc * ( NX * deg + 1 )
338 ofstream ltgout(
"ltg.dat");
339 for (
int k=0;k<NZ;k++)
341 for (
int j=0; j<NY; j++)
343 for (
int i=0; i<NX; i++)
345 ielem = i + j * NX + k * NX * NY;
346 for (
int m=0; m<numLineFieldsG*numLineFieldsG*numLineFieldsG; m++)
348 ltgout << ltgMapping(ielem,m) <<
" ";
358 Epetra_SerialComm Comm;
359 Epetra_Map globalMapG(numDOF, 0, Comm);
360 Epetra_FEVector u(globalMapG); u.Random();
361 Epetra_FEVector Ku(globalMapG);
370 Epetra_Time multTimer(Comm);
371 Teuchos::BLAS<int,double> blas;
375 double *uVals = u[0];
376 double *KuVals = Ku[0];
378 Epetra_Time scatterTimer(Comm);
379 std::cout <<
"Scattering\n";
381 for (
int k=0; k<numElems; k++)
383 for (
int i=0;i<numLineFieldsG*numLineFieldsG*numLineFieldsG;i++)
385 uScattered(k,i) = uVals[ltgMapping(k,i)];
390 const double scatterTime = scatterTimer.ElapsedTime();
396 Epetra_Time localAppTimer(Comm);
398 uScattered.resize(numElems,numLineFieldsG,numLineFieldsG,numLineFieldsG);
404 for (ielem=0;ielem<numElems;ielem++)
409 for (
int k=0;k<numLineFieldsG;k++)
411 for (
int j=0;j<numLineFieldsG;j++)
413 for (
int i=0;i<numLineFieldsG;i++)
423 for (
int k=0;k<numLineFieldsG;k++)
425 for (
int j=0;j<numLineFieldsG;j++)
427 for (
int i=0;i<numLineFieldsG;i++)
429 for (
int q=0;q<numLineFieldsG;q++)
431 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(i,q,0);
440 for (
int k=0;k<numLineFieldsG;k++)
442 const double wt1 = hcur * cubWeights1D(k);
443 for (
int j=0;j<numLineFieldsG;j++)
445 const double wtstuff = wt1 * cubWeights1D(j);
446 for (
int i=0;i<numLineFieldsG;i++)
448 for (
int q=0;q<numLineFieldsG;q++)
450 KuScattered(ielem,cur) += wtstuff
451 * cubWeights1D(q) * Du(k,j,q) * lineGrads(i,q,0);
462 for (
int k=0;k<numLineFieldsG;k++)
464 for (
int j=0;j<numLineFieldsG;j++)
466 for (
int i=0;i<numLineFieldsG;i++)
474 for (
int k=0;k<numLineFieldsG;k++)
476 for (
int j=0;j<numLineFieldsG;j++)
478 for (
int i=0;i<numLineFieldsG;i++)
480 for (
int q=0;q<numLineFieldsG;q++)
482 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(j,q,0);
491 for (
int k=0;k<numLineFieldsG;k++)
493 const double wt1 = hcur * cubWeights1D(k);
494 for (
int j=0;j<numLineFieldsG;j++)
496 for (
int i=0;i<numLineFieldsG;i++)
498 const double wtstuff = cubWeights1D(i) * wt1;
499 for (
int q=0;q<numLineFieldsG;q++)
501 KuScattered(ielem,cur) += wtstuff * cubWeights1D(q) * Du(k,q,i) * lineGrads(j,q,0);
512 for (
int k=0;k<numLineFieldsG;k++)
514 for (
int j=0;j<numLineFieldsG;j++)
516 for (
int i=0;i<numLineFieldsG;i++)
524 for (
int k=0;k<numLineFieldsG;k++)
526 for (
int j=0;j<numLineFieldsG;j++)
528 for (
int i=0;i<numLineFieldsG;i++)
530 for (
int q=0;q<numLineFieldsG;q++)
532 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(k,q,0);
541 for (
int k=0;k<numLineFieldsG;k++)
543 for (
int j=0;j<numLineFieldsG;j++)
545 const double wt1 = hcur * cubWeights1D(j);
546 for (
int i=0;i<numLineFieldsG;i++)
548 const double wtstuff = cubWeights1D(i) * wt1;
549 for (
int q=0;q<numLineFieldsG;q++)
551 KuScattered(ielem,cur) += wtstuff * cubWeights1D(q) * Du(q,j,i) * lineGrads(k,q,0);
560 const double localAppTime = localAppTimer.ElapsedTime();
562 Epetra_Time gatherTimer(Comm);
564 for (
int k=0;k<numElems;k++)
566 for (
int i=0;i<numLineFieldsG*numLineFieldsG*numLineFieldsG;i++)
568 KuVals[ltgMapping(k,i)] += KuScattered(k,i);
572 const double gatherTime = gatherTimer.ElapsedTime();
575 *outStream <<
"Time to scatter " << scatterTime <<
"\n";
576 *outStream <<
"Time for local application " << localAppTime <<
"\n";
577 *outStream <<
"Time to gather " << gatherTime <<
"\n";
578 *outStream <<
"Total matrix-free time " << scatterTime + localAppTime + gatherTime <<
"\n";
581 *outStream <<
"End Result: TEST PASSED\n";
584 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubaturePolylib class.
Header file for utility class to provide multidimensional containers.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....