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
85 #include "Intrepid_CellTools.hpp"
86 #include "Intrepid_ArrayTools.hpp"
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 
111 using namespace std;
112 using namespace Intrepid;
113 
114 int 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  int 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  int 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 
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid utilities.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
A factory class that generates specific instances of cubatures.
Header file for the Intrepid::CubaturePolylib class.