rect2DMeshGenerator.cpp

Go to the documentation of this file.
00001 #include "rect2DMeshGenerator.hpp"
00002 #include "Epetra_IntSerialDenseVector.h"
00003 #include "Epetra_SerialDenseMatrix.h"
00004 #include "Epetra_IntSerialDenseMatrix.h"
00005 #include "Teuchos_TestForException.hpp"
00006 #include "Teuchos_FancyOStream.hpp"
00007 
00008 void GLpApp::rect2DMeshGenerator(
00009   const int                      numProc
00010   ,const int                     procRank
00011   ,const double                  len_x
00012   ,const double                  len_y
00013   ,const int                     local_nx
00014   ,const int                     local_ny
00015   ,const int                     bndy_marker
00016   ,Epetra_IntSerialDenseVector   *ipindx_out
00017   ,Epetra_SerialDenseMatrix      *ipcoords_out
00018   ,Epetra_IntSerialDenseVector   *pindx_out
00019   ,Epetra_SerialDenseMatrix      *pcoords_out
00020   ,Epetra_IntSerialDenseMatrix   *t_out
00021   ,Epetra_IntSerialDenseMatrix   *e_out
00022   ,std::ostream                  *out_arg
00023   ,const bool                    dumpAll
00024   )
00025 {
00026   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00027     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00028   Teuchos::OSTab tab(out);
00029   if(out.get())
00030     *out << "\nGenerating rectangular mesh with len_x="<<len_x<<", len_y="<<len_y<<", local_nx="<<local_nx<<", and local_ny="<<local_ny<<" ...\n";
00031   //
00032   // Validate input
00033   //
00034 #ifdef TEUCHOS_DEBUG
00035   TEST_FOR_EXCEPT(len_x <= 0.0);
00036   TEST_FOR_EXCEPT(len_y <= 0.0);
00037   TEST_FOR_EXCEPT(local_nx <= 0);
00038   TEST_FOR_EXCEPT(local_ny <= 0);
00039   TEST_FOR_EXCEPT(ipindx_out==NULL);
00040   TEST_FOR_EXCEPT(ipcoords_out==NULL);
00041   TEST_FOR_EXCEPT(pindx_out==NULL);
00042   TEST_FOR_EXCEPT(pcoords_out==NULL);
00043   TEST_FOR_EXCEPT(t_out==NULL);
00044   TEST_FOR_EXCEPT(e_out==NULL);
00045 #endif
00046   //
00047   // Get local references
00048   //
00049   Epetra_IntSerialDenseVector   &ipindx = *ipindx_out;
00050   Epetra_SerialDenseMatrix      &ipcoords = *ipcoords_out;
00051   Epetra_IntSerialDenseVector   &pindx = *pindx_out;
00052   Epetra_SerialDenseMatrix      &pcoords = *pcoords_out;
00053   Epetra_IntSerialDenseMatrix   &t = *t_out;
00054   Epetra_IntSerialDenseMatrix   &e = *e_out;
00055   //
00056   // Get dimensions
00057   //
00058   const int
00059     numip = ( local_nx + 1 ) * ( local_ny + ( procRank == numProc-1 ? 1 : 0 ) ),
00060     numcp = ( procRank == numProc-1 ? 0 : (local_nx+1) ),
00061     nump  = numip + numcp,
00062     numelems = 2*local_nx*local_ny,
00063     numedges = 2*local_ny + ( procRank==0 ? local_nx : 0 ) + ( procRank==numProc-1 ? local_nx : 0 );
00064   const double
00065     delta_x = len_x / local_nx,
00066     delta_y = (len_y/numProc) / local_ny;
00067   if(out.get()) {
00068     *out
00069       << "\nnumip = " << numip
00070       << "\nnumcp = " << numcp
00071       << "\nnump = " << nump
00072       << "\nnumelems = " << numelems
00073       << "\nnumedges = " << numedges
00074       << "\ndelta_x = " << delta_x
00075       << "\ndelta_y = " << delta_y
00076       << "\n";
00077   }
00078   //
00079   // Resize mesh storage
00080   //
00081   ipindx.Size(numip);
00082   ipcoords.Shape(numip,2);
00083   pindx.Size(nump);
00084   pcoords.Shape(nump,2);
00085   t.Shape(numelems,3);
00086   e.Shape(numedges,3);
00087   //
00088   // Get the global offsets for the local nodes IDs and y coordinates
00089   //
00090   const int     localNodeOffset   = procRank*(local_nx+1)*local_ny;
00091   const double  localYCoordOffset = procRank*delta_y*local_ny;
00092   //
00093   // Set the local node IDs and their cooridinates
00094   //
00095   if(out.get()) *out << "\nGenerating the node list ...\n";
00096   tab.incrTab();
00097   // Owned and shared
00098   int node_i  = 0;
00099   int global_node_id = localNodeOffset+1;
00100   for( int i_y = 0; i_y < local_ny + 1; ++i_y ) {
00101     for( int i_x = 0; i_x < local_nx + 1; ++i_x ) {
00102       pindx(node_i) = global_node_id;
00103       pcoords(node_i,0) = i_x * delta_x;
00104       pcoords(node_i,1) = i_y * delta_y + localYCoordOffset;
00105       if(out.get() && dumpAll)
00106         *out << "node_i="<<node_i<<",global_node_id="<<global_node_id<<",("<<pcoords(node_i,0)<<","<<pcoords(node_i,1)<<")\n";
00107       ++node_i;
00108       ++global_node_id;
00109     }
00110   }
00111   tab.incrTab(-1);
00112   TEST_FOR_EXCEPT(node_i != nump);
00113   // Locally owned only
00114   for( int i = 0; i < numip; ++i ) {
00115     ipindx(i) = pindx(i);
00116     ipcoords(i,0) = pcoords(i,0);
00117     ipcoords(i,1) = pcoords(i,1);
00118   }
00119   //
00120   // Set the elements
00121   //
00122   if(out.get()) *out << "\nGenerating the element list ...\n";
00123   tab.incrTab();
00124   int ele_k = 0;
00125   global_node_id = localNodeOffset+1;
00126   for( int i_y = 0; i_y < local_ny; ++i_y ) {
00127     for( int i_x = 0; i_x < local_nx; ++i_x ) {
00128       t(ele_k,0) = global_node_id;
00129       t(ele_k,1) = global_node_id + (local_nx + 1);
00130       t(ele_k,2) = global_node_id + 1;
00131       if(out.get() && dumpAll)
00132         *out << "ele_k="<<ele_k<<",("<<t(ele_k,0)<<","<<t(ele_k,1)<<","<<t(ele_k,2)<<")\n";
00133       ++ele_k;
00134       t(ele_k,0) = global_node_id + 1;
00135       t(ele_k,1) = global_node_id + (local_nx + 1);
00136       t(ele_k,2) = global_node_id + (local_nx + 1) + 1;
00137       if(out.get() && dumpAll)
00138         *out << "ele_k="<<ele_k<<",("<<t(ele_k,0)<<","<<t(ele_k,1)<<","<<t(ele_k,2)<<")\n";
00139       ++ele_k;
00140       ++global_node_id;
00141     }
00142     ++global_node_id;
00143   }
00144   tab.incrTab(-1);
00145   TEST_FOR_EXCEPT(ele_k != numelems);
00146   //
00147   // Set the edges
00148   //
00149   int edge_j = 0;
00150   if(procRank==0) {
00151     // Bottom edges
00152     if(out.get()) *out << "\nGenerating the bottom edges ...\n";
00153     tab.incrTab();
00154     global_node_id = localNodeOffset+1;
00155     for( int i_x = 0; i_x < local_nx; ++i_x ) {
00156       e(edge_j,0) = global_node_id;
00157       e(edge_j,1) = global_node_id + 1;
00158       e(edge_j,2) = bndy_marker;
00159       if(out.get() && dumpAll)
00160         *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
00161       ++edge_j;
00162       global_node_id += 1;
00163     }
00164     tab.incrTab(-1);
00165   }
00166   // Left edges
00167   if(out.get()) *out << "\nGenerating the left edges ...\n";
00168   tab.incrTab();
00169   global_node_id = localNodeOffset+1;
00170   for( int i_y = 0; i_y < local_ny; ++i_y ) {
00171     e(edge_j,0) = global_node_id;
00172     e(edge_j,1) = global_node_id + (local_nx + 1);
00173     e(edge_j,2) = bndy_marker;
00174     if(out.get() && dumpAll)
00175       *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
00176     ++edge_j;
00177     global_node_id += (local_nx + 1);
00178   }
00179   tab.incrTab(-1);
00180   // Right edges
00181   if(out.get()) *out << "\nGenerating the right edges ...\n";
00182   tab.incrTab();
00183   global_node_id = localNodeOffset + 1 + local_nx;
00184   for( int i_y = 0; i_y < local_ny; ++i_y ) {
00185     e(edge_j,0) = global_node_id;
00186     e(edge_j,1) = global_node_id + (local_nx + 1);
00187     e(edge_j,2) = bndy_marker;
00188     if(out.get() && dumpAll)
00189       *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
00190     ++edge_j;
00191     global_node_id += (local_nx + 1);
00192   }
00193   tab.incrTab(-1);
00194   if(procRank==numProc-1) {
00195     // Top edges
00196     if(out.get()) *out << "\nGenerating the top edges ...\n";
00197     tab.incrTab();
00198     global_node_id = localNodeOffset+1+(local_nx+1)*local_ny;
00199     for( int i_x = 0; i_x < local_nx; ++i_x ) {
00200       e(edge_j,0) = global_node_id;
00201       e(edge_j,1) = global_node_id + 1;
00202       e(edge_j,2) = bndy_marker;
00203       if(out.get() && dumpAll)
00204         *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
00205       ++edge_j;
00206       global_node_id += 1;
00207     }
00208     tab.incrTab(-1);
00209   }
00210   TEST_FOR_EXCEPT(edge_j != numedges);
00211 }

Generated on Thu Sep 18 12:31:44 2008 for EpetraExt by doxygen 1.3.9.1