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