EpetraExt Development
Loading...
Searching...
No Matches
rect2DMeshGenerator.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services Package
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#include "rect2DMeshGenerator.hpp"
45#include "Epetra_IntSerialDenseVector.h"
46#include "Epetra_SerialDenseMatrix.h"
47#include "Epetra_IntSerialDenseMatrix.h"
48#include "Teuchos_Assert.hpp"
49#include "Teuchos_FancyOStream.hpp"
50#include "Teuchos_RCP.hpp"
51
52void GLpApp::rect2DMeshGenerator(
53 const int numProc
54 ,const int procRank
55 ,const double len_x
56 ,const double len_y
57 ,const int local_nx
58 ,const int local_ny
59 ,const int bndy_marker
60 ,Epetra_IntSerialDenseVector *ipindx_out
61 ,Epetra_SerialDenseMatrix *ipcoords_out
62 ,Epetra_IntSerialDenseVector *pindx_out
63 ,Epetra_SerialDenseMatrix *pcoords_out
64 ,Epetra_IntSerialDenseMatrix *t_out
65 ,Epetra_IntSerialDenseMatrix *e_out
66 ,std::ostream *out_arg
67 ,const bool dumpAll
68 )
69{
70 Teuchos::RCP<Teuchos::FancyOStream>
71 out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
72 Teuchos::OSTab tab(out);
73 if(out.get())
74 *out << "\nGenerating rectangular mesh with len_x="<<len_x<<", len_y="<<len_y<<", local_nx="<<local_nx<<", and local_ny="<<local_ny<<" ...\n";
75 //
76 // Validate input
77 //
78#ifdef TEUCHOS_DEBUG
79 TEUCHOS_TEST_FOR_EXCEPT(len_x <= 0.0);
80 TEUCHOS_TEST_FOR_EXCEPT(len_y <= 0.0);
81 TEUCHOS_TEST_FOR_EXCEPT(local_nx <= 0);
82 TEUCHOS_TEST_FOR_EXCEPT(local_ny <= 0);
83 TEUCHOS_TEST_FOR_EXCEPT(ipindx_out==NULL);
84 TEUCHOS_TEST_FOR_EXCEPT(ipcoords_out==NULL);
85 TEUCHOS_TEST_FOR_EXCEPT(pindx_out==NULL);
86 TEUCHOS_TEST_FOR_EXCEPT(pcoords_out==NULL);
87 TEUCHOS_TEST_FOR_EXCEPT(t_out==NULL);
88 TEUCHOS_TEST_FOR_EXCEPT(e_out==NULL);
89#endif
90 //
91 // Get local references
92 //
93 Epetra_IntSerialDenseVector &ipindx = *ipindx_out;
94 Epetra_SerialDenseMatrix &ipcoords = *ipcoords_out;
95 Epetra_IntSerialDenseVector &pindx = *pindx_out;
96 Epetra_SerialDenseMatrix &pcoords = *pcoords_out;
97 Epetra_IntSerialDenseMatrix &t = *t_out;
98 Epetra_IntSerialDenseMatrix &e = *e_out;
99 //
100 // Get dimensions
101 //
102 const int
103 numip = ( local_nx + 1 ) * ( local_ny + ( procRank == numProc-1 ? 1 : 0 ) ),
104 numcp = ( procRank == numProc-1 ? 0 : (local_nx+1) ),
105 nump = numip + numcp,
106 numelems = 2*local_nx*local_ny,
107 numedges = 2*local_ny + ( procRank==0 ? local_nx : 0 ) + ( procRank==numProc-1 ? local_nx : 0 );
108 const double
109 delta_x = len_x / local_nx,
110 delta_y = (len_y/numProc) / local_ny;
111 if(out.get()) {
112 *out
113 << "\nnumip = " << numip
114 << "\nnumcp = " << numcp
115 << "\nnump = " << nump
116 << "\nnumelems = " << numelems
117 << "\nnumedges = " << numedges
118 << "\ndelta_x = " << delta_x
119 << "\ndelta_y = " << delta_y
120 << "\n";
121 }
122 //
123 // Resize mesh storage
124 //
125 ipindx.Size(numip);
126 ipcoords.Shape(numip,2);
127 pindx.Size(nump);
128 pcoords.Shape(nump,2);
129 t.Shape(numelems,3);
130 e.Shape(numedges,3);
131 //
132 // Get the global offsets for the local nodes IDs and y coordinates
133 //
134 const int localNodeOffset = procRank*(local_nx+1)*local_ny;
135 const double localYCoordOffset = procRank*delta_y*local_ny;
136 //
137 // Set the local node IDs and their cooridinates
138 //
139 if(out.get()) *out << "\nGenerating the node list ...\n";
140 tab.incrTab();
141 // Owned and shared
142 int node_i = 0;
143 int global_node_id = localNodeOffset+1;
144 for( int i_y = 0; i_y < local_ny + 1; ++i_y ) {
145 for( int i_x = 0; i_x < local_nx + 1; ++i_x ) {
146 pindx(node_i) = global_node_id;
147 pcoords(node_i,0) = i_x * delta_x;
148 pcoords(node_i,1) = i_y * delta_y + localYCoordOffset;
149 if(out.get() && dumpAll)
150 *out << "node_i="<<node_i<<",global_node_id="<<global_node_id<<",("<<pcoords(node_i,0)<<","<<pcoords(node_i,1)<<")\n";
151 ++node_i;
152 ++global_node_id;
153 }
154 }
155 tab.incrTab(-1);
156 TEUCHOS_TEST_FOR_EXCEPT(node_i != nump);
157 // Locally owned only
158 for( int i = 0; i < numip; ++i ) {
159 ipindx(i) = pindx(i);
160 ipcoords(i,0) = pcoords(i,0);
161 ipcoords(i,1) = pcoords(i,1);
162 }
163 //
164 // Set the elements
165 //
166 if(out.get()) *out << "\nGenerating the element list ...\n";
167 tab.incrTab();
168 int ele_k = 0;
169 global_node_id = localNodeOffset+1;
170 for( int i_y = 0; i_y < local_ny; ++i_y ) {
171 for( int i_x = 0; i_x < local_nx; ++i_x ) {
172 t(ele_k,0) = global_node_id;
173 t(ele_k,1) = global_node_id + (local_nx + 1);
174 t(ele_k,2) = global_node_id + 1;
175 if(out.get() && dumpAll)
176 *out << "ele_k="<<ele_k<<",("<<t(ele_k,0)<<","<<t(ele_k,1)<<","<<t(ele_k,2)<<")\n";
177 ++ele_k;
178 t(ele_k,0) = global_node_id + 1;
179 t(ele_k,1) = global_node_id + (local_nx + 1);
180 t(ele_k,2) = global_node_id + (local_nx + 1) + 1;
181 if(out.get() && dumpAll)
182 *out << "ele_k="<<ele_k<<",("<<t(ele_k,0)<<","<<t(ele_k,1)<<","<<t(ele_k,2)<<")\n";
183 ++ele_k;
184 ++global_node_id;
185 }
186 ++global_node_id;
187 }
188 tab.incrTab(-1);
189 TEUCHOS_TEST_FOR_EXCEPT(ele_k != numelems);
190 //
191 // Set the edges
192 //
193 int edge_j = 0;
194 if(procRank==0) {
195 // Bottom edges
196 if(out.get()) *out << "\nGenerating the bottom edges ...\n";
197 tab.incrTab();
198 global_node_id = localNodeOffset+1;
199 for( int i_x = 0; i_x < local_nx; ++i_x ) {
200 e(edge_j,0) = global_node_id;
201 e(edge_j,1) = global_node_id + 1;
202 e(edge_j,2) = bndy_marker;
203 if(out.get() && dumpAll)
204 *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
205 ++edge_j;
206 global_node_id += 1;
207 }
208 tab.incrTab(-1);
209 }
210 // Left edges
211 if(out.get()) *out << "\nGenerating the left edges ...\n";
212 tab.incrTab();
213 global_node_id = localNodeOffset+1;
214 for( int i_y = 0; i_y < local_ny; ++i_y ) {
215 e(edge_j,0) = global_node_id;
216 e(edge_j,1) = global_node_id + (local_nx + 1);
217 e(edge_j,2) = bndy_marker;
218 if(out.get() && dumpAll)
219 *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
220 ++edge_j;
221 global_node_id += (local_nx + 1);
222 }
223 tab.incrTab(-1);
224 // Right edges
225 if(out.get()) *out << "\nGenerating the right edges ...\n";
226 tab.incrTab();
227 global_node_id = localNodeOffset + 1 + local_nx;
228 for( int i_y = 0; i_y < local_ny; ++i_y ) {
229 e(edge_j,0) = global_node_id;
230 e(edge_j,1) = global_node_id + (local_nx + 1);
231 e(edge_j,2) = bndy_marker;
232 if(out.get() && dumpAll)
233 *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
234 ++edge_j;
235 global_node_id += (local_nx + 1);
236 }
237 tab.incrTab(-1);
238 if(procRank==numProc-1) {
239 // Top edges
240 if(out.get()) *out << "\nGenerating the top edges ...\n";
241 tab.incrTab();
242 global_node_id = localNodeOffset+1+(local_nx+1)*local_ny;
243 for( int i_x = 0; i_x < local_nx; ++i_x ) {
244 e(edge_j,0) = global_node_id;
245 e(edge_j,1) = global_node_id + 1;
246 e(edge_j,2) = bndy_marker;
247 if(out.get() && dumpAll)
248 *out << "edge_j="<<edge_j<<",("<<e(edge_j,0)<<","<<e(edge_j,1)<<"),"<<e(edge_j,2)<<"\n";
249 ++edge_j;
250 global_node_id += 1;
251 }
252 tab.incrTab(-1);
253 }
254 TEUCHOS_TEST_FOR_EXCEPT(edge_j != numedges);
255}
Epetra_IntSerialDenseMatrix
Epetra_IntSerialDenseMatrix::Shape
int Shape(int NumRows, int NumCols)
Epetra_IntSerialDenseVector
Epetra_IntSerialDenseVector::Size
int Size(int Length_in)
Epetra_SerialDenseMatrix
Epetra_SerialDenseMatrix::Shape
int Shape(int NumRows, int NumCols)
GLpApp::rect2DMeshGenerator
void rect2DMeshGenerator(const int numProc, const int procRank, const double len_x, const double len_y, const int local_nx, const int local_ny, const int bndy_marker, Epetra_IntSerialDenseVector *ipindx_out, Epetra_SerialDenseMatrix *ipcoords_out, Epetra_IntSerialDenseVector *pindx_out, Epetra_SerialDenseMatrix *pcoords_out, Epetra_IntSerialDenseMatrix *t_out, Epetra_IntSerialDenseMatrix *e_out, std::ostream *out, const bool dumpAll)
Generate a simple rectangular 2D triangular mesh that is only partitioned between processors in the y...
Definition rect2DMeshGenerator.cpp:52
Generated by doxygen 1.9.8