NeoPZ
tpzhierarquicalgrid.cpp
Go to the documentation of this file.
1 /*
2  <one line to give the program's name and a brief idea of what it does.>
3  Copyright (C) 2014 <copyright holder> <email>
4 
5  This program is free software: you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 
20 #include "tpzhierarquicalgrid.h"
21 #include "pzgmesh.h"
22 #include "pzgeoel.h"
23 #include "TPZGeoLinear.h"
24 #include "pzgeopoint.h"
25 #include "tpztriangle.h"
26 #include "pzgeoquad.h"
27 #include "TPZRefPattern.h"
28 #include "tpzgeoelrefpattern.h"
29 #include "TPZGeoCube.h"
30 #include "pzgeoprism.h"
31 #include "pzgeotetrahedra.h"
32 
34 {
35 
36  fFileName = "untitled";
37  fNonAffineQ = false;
38  fIsQuad = true;
39  fIsPrism = false;
40  fIsTetrahedron = false;
41  fComputedGeomesh = NULL;
42  ffrontMatID = -999;
43  fbackMatID = -1000;
44  fSubBases.Resize(0);
45  fBase = NULL;
46 }
47 
49 {
50  if(Geomesh->NElements() == 0)
51  {
52  if(fBase) std::cout << "Number of elements" << fBase->NElements() << std::endl;
53  DebugStop();
54  }
55 
56  if(Geomesh->NNodes() == 0)
57  {
58  std::cout << "Number of nodes" << fBase->NElements() << std::endl;
59  DebugStop();
60  }
61 
62  fFileName = "untitled";
63  fNonAffineQ = false;
64  fIsQuad = true;
65  fIsPrism = false;
66  fIsTetrahedron = false;
67  fComputedGeomesh = NULL;
68  ffrontMatID = -999;
69  fbackMatID = -1000;
70  fBase = Geomesh;
71  fSubBases.Resize(0);
72 }
73 
74 
75 
77 {
78 
79  fFileName = other.fFileName;
80  fNonAffineQ = other.fNonAffineQ;
81  fIsQuad = other.fIsQuad;
82  fIsPrism = other.fIsPrism;
85  ffrontMatID = other.ffrontMatID;
86  fbackMatID = other.fbackMatID;
87  fSubBases = other.fSubBases;
88  fBase = other.fBase;
89 
90 }
91 
93 {
94 
95 }
96 
98 {
99  fFileName = other.fFileName;
100  fNonAffineQ = other.fNonAffineQ;
101  fIsQuad = other.fIsQuad;
102  fIsPrism = other.fIsPrism;
105  ffrontMatID = other.ffrontMatID;
106  fbackMatID = other.fbackMatID;
107  fSubBases = other.fSubBases;
108  fBase = other.fBase;
109  return *this;
110 }
111 
113 {
114  DebugStop();
116  return 0; //to fix WIN compiler error
117 }
118 
120 {
121  fSubBases.Resize(n+1);
123  int NNodesBase = fBase->NNodes();
124  fComputedGeomesh->NodeVec().Resize( (n+1) * NNodesBase);
125 
126  // Creating new elements
127  int nodeId = 0;
128  REAL sing = 1.0;
129  for(int il = 0; il < (n+1); il++ )
130  {
131  // copying l extrusions
132  // For a while all of them are equal
133  // fSubBases[il] = new TPZGeoMesh(fBase);
134 
135  for(int inode = 0; inode < NNodesBase; inode++)
136  {
137  TPZVec<REAL> tpara(1),NewCoordinates_r(3,0.0);
138  TPZVec<REAL> NewCoordinates(3,0.0);
139  tpara[0] = dt*il+t;
140 
141  TPZVec<REAL> Coordinates(3,0.0);
142  fComputedGeomesh->NodeVec()[inode + il * NNodesBase] = fBase->NodeVec()[inode];
143  // NewGeomesh->NodeVec()[inode + il*fBase->NNodes()] = fSubBases[il]->NodeVec()[inode];
144  fComputedGeomesh->NodeVec()[inode + il * NNodesBase].GetCoordinates(Coordinates);
145 
146  fParametricFunction->Execute(tpara,NewCoordinates);
147  NewCoordinates_r[0] = REAL(NewCoordinates[0]);
148  NewCoordinates_r[1] = REAL(NewCoordinates[1]);
149  NewCoordinates_r[2] = REAL(NewCoordinates[2]);
150 
151  if(((il+1)%2==0 && fNonAffineQ) && fBase->Dimension() == 2){
152  Coordinates[0]+=NewCoordinates_r[0];
153  Coordinates[1]+=NewCoordinates_r[1];
154  Coordinates[2]+=NewCoordinates_r[2];
155  Coordinates[2]+= sing*dt/2.0;
156  sing *= -1.0;
157  }
158  else{
159  Coordinates[0]+=NewCoordinates_r[0];
160  Coordinates[1]+=NewCoordinates_r[1];
161  Coordinates[2]+=NewCoordinates_r[2];
162  }
163 
164  fComputedGeomesh->NodeVec()[inode + il * NNodesBase].SetCoord(Coordinates);
165  fComputedGeomesh->NodeVec()[inode + il * NNodesBase].SetNodeId(nodeId);
166  nodeId++;
167  }
168  }
169 
171 
172 // int fbasedim = fBase->Dimension();
173 
174  int elid=0;
175  for(int iel = 0; iel < fBase->NElements(); iel++)
176  {
177  int ielDim = fBase->ElementVec()[iel]->Dimension();
178  int ielMatId = fBase->ElementVec()[iel]->MaterialId();
179  CreateGeometricElement(n,iel,ielDim,ielMatId,elid);
180  }
181 
184  return fComputedGeomesh;
185 
186 }
187 
188 void TPZHierarquicalGrid::CreateGeometricElement(int n, int iel,int eldim, int elmatid, int &elid)
189 {
190  int jump = fBase->NNodes();
191  int dim = fBase->Dimension();
192 
193  TPZGeoEl *gel = fBase->ElementVec()[iel];
194  int gelNodes = gel->NNodes();
195 
196  // Computing current topology
197  TPZManVector<int64_t,10> CTopology(gelNodes);
198  for(int inode = 0; inode < CTopology.size(); inode++)
199  {
200  CTopology[inode] = gel->Node(inode).Id();
201  }
202 
203  bool Is1D = false;
204  bool Is2D = false;
205  bool Is3D = false;
206 
207  for(int il = 1; il < (n+1); il++ )
208  {
209 
210  switch (eldim) {
211  case 0:
212  {
213 
214  //Defining boundaries
215  if (dim==0) {
216  if (il==1){
217  TPZVec<int64_t> Topology(gelNodes);
218  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
220  }
221  if (il==n)
222  {
223  TPZVec<int64_t> Topology(gelNodes);
224  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
226  }
227  }
228 
229  // if (dim==1) {
230  // if (il==1){
231  // TPZVec<int64_t> Topology(gelNodes+1);
232  // Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
233  // Topology[1]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
234  // new TPZGeoElRefPattern < pzgeom::TPZGeoLinear > (elid++, Topology, elmatid,*fComputedGeomesh);
235  // }
236  // if (il==n)
237  // {
238  // TPZVec<int64_t> Topology(gelNodes+1);
239  // Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
240  // Topology[1]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
241  // new TPZGeoElRefPattern < pzgeom::TPZGeoLinear > (elid++, Topology, elmatid,*fComputedGeomesh);
242  //
243  // }
244  // }
245 
246  TPZVec<int64_t> Topology(gelNodes+1);
247  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
248  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
249  new TPZGeoElRefPattern < pzgeom::TPZGeoLinear > (elid++, Topology, elmatid,*fComputedGeomesh);
250  Is1D = true;
251  }
252  break;
253  case 1:
254  {
255  if (dim==1) {
256  if (il==1){
257  TPZVec<int64_t> Topology(gelNodes+1);
258  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
259  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
261  }
262  if (il==n)
263  {
264  TPZVec<int64_t> Topology(gelNodes+1);
265  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
266  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
268  }
269  }
270 
271  if (fIsQuad) {
272  // quadrilateras
273  TPZVec<int64_t> Topology(gelNodes+2);
274  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
275  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
276  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
277  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
278  new TPZGeoElRefPattern < pzgeom::TPZGeoQuad > (elid++, Topology, elmatid,*fComputedGeomesh);
279  }
280  else
281  {
282  // triangles
283  TPZVec<int64_t> Topology(gelNodes+1);
284  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
285  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
286  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
287  new TPZGeoElRefPattern < pzgeom::TPZGeoTriangle > (elid++, Topology, elmatid,*fComputedGeomesh);
288 
289  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
290  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
291  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
292  new TPZGeoElRefPattern < pzgeom::TPZGeoTriangle > (elid++, Topology, elmatid,*fComputedGeomesh);
293 
294  }
295  Is2D = true;
296 
297  }
298  break;
299  case 2:
300  {
301 
302  if (dim==2) {
303  if (il==1){
304 
305  if (gel->Type() == EQuadrilateral) {
306  // quadrilateras
307  TPZVec<int64_t> Topology(gelNodes+2);
308  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
309  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
310  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
311  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[3] + (il - 1) * jump].Id();
313  }
314  else{
315  // triangles
316  TPZVec<int64_t> Topology(gelNodes+1);
317  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
318  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
319  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
321 
322  }
323  }
324  if (il==n)
325  {
326  if (gel->Type() == EQuadrilateral) {
327  // quadrilateras
328  TPZVec<int64_t> Topology(gelNodes+2);
329  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
330  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
331  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
332  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[3] + (il - 0) * jump].Id();
334  }
335  else{
336  // triangles
337  TPZVec<int64_t> Topology(gelNodes+1);
338  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
339  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
340  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
342 
343  }
344  }
345  }
346 
347  if (!fIsTetrahedron) {
348 
349  if (fIsPrism) {
350  // Prisms
351  TPZVec<int64_t> Topology(2*gelNodes);
352  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
353  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
354  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
355  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
356  Topology[4]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
357  Topology[5]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
358  new TPZGeoElRefPattern < pzgeom::TPZGeoPrism > (elid++, Topology, elmatid,*fComputedGeomesh);
359  }
360  else
361  {
362  // Cubes
363  TPZVec<int64_t> Topology(gelNodes+4);
364  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
365  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
366  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
367  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[3] + (il - 1) * jump].Id();
368  Topology[4]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
369  Topology[5]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
370  Topology[6]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
371  Topology[7]=fComputedGeomesh->NodeVec()[CTopology[3] + (il - 0) * jump].Id();
372  new TPZGeoElRefPattern < pzgeom::TPZGeoCube > (elid++, Topology, elmatid,*fComputedGeomesh);
373  }
374  }
375  else
376  {
377  if (fIsPrism) {
378  // Prisms
379  TPZVec<int64_t> Topology(2*gelNodes);
380  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
381  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
382  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
383  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
384  Topology[4]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
385  Topology[5]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
386  new TPZGeoElRefPattern < pzgeom::TPZGeoPrism > (elid++, Topology, elmatid,*fComputedGeomesh);
387  }
388  else{
389  // Tetrahedron
390  TPZVec<int64_t> Topology(gelNodes+1);
391  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 1) * jump].Id();
392  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
393  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
394  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
395  new TPZGeoElRefPattern < pzgeom::TPZGeoTetrahedra > (elid++, Topology, elmatid,*fComputedGeomesh);
396 
397  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 1) * jump].Id();
398  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
399  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
400  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
401  new TPZGeoElRefPattern < pzgeom::TPZGeoTetrahedra > (elid++, Topology, elmatid,*fComputedGeomesh);
402 
403  Topology[0]=fComputedGeomesh->NodeVec()[CTopology[0] + (il - 0) * jump].Id();
404  Topology[1]=fComputedGeomesh->NodeVec()[CTopology[1] + (il - 0) * jump].Id();
405  Topology[2]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 0) * jump].Id();
406  Topology[3]=fComputedGeomesh->NodeVec()[CTopology[2] + (il - 1) * jump].Id();
407  new TPZGeoElRefPattern < pzgeom::TPZGeoTetrahedra > (elid++, Topology, elmatid,*fComputedGeomesh);
408 
409  }
410 
411  }
412  Is3D = true;
413 
414  }
415  break;
416  default:
417  {
418  std::cout << "Connection not implemented " << std::endl;
419  DebugStop();
420  }
421  break;
422  }
423  }
424 
425  if (Is1D) {
427  }
428 
429  if (Is2D) {
431  }
432 
433  if (Is3D) {
435  }
436 
437 }
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
void SetMaxElementId(int64_t id)
Used in patch meshes.
Definition: pzgmesh.h:123
bool operator==(const TPZHierarquicalGrid &other) const
TPZAutoPointer< TPZGeoMesh > fBase
A geometric mesh generated from other sources.
TPZAutoPointer< TPZFunction< REAL > > fParametricFunction
Pointer to parametric function of t parameter.
int fbackMatID
Extrusion back material id.
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
bool fIsQuad
2d extrusion is based on quadrilaterals
int64_t NElements() const
Number of elements of the mesh.
Definition: pzgmesh.h:129
std::string fFileName
Name of the fine mesh to be extended.
void CreateGeometricElement(int n, int iel, int eldim, int elmatid, int &elid)
void SetDimension(int dim)
Set Dimension.
Definition: pzgmesh.h:285
Contains the TPZTriangle class which defines the topology of a triangle.
TPZGeoMesh * ComputeExtrusion(REAL t, REAL dt, int n)
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
void Resize(const int newsize)
Increase the size of the chunk vector.
Definition: pzadmchunk.h:280
Contains declaration of TPZMesh class which defines a geometrical mesh and contains a corresponding l...
void SetMaxNodeId(int64_t id)
Used in patch meshes.
Definition: pzgmesh.h:120
Implements a generic geometric element which is refined according to a generic refinement pattern...
Definition: pzgmesh.h:35
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
Contains the TPZGeoTetrahedra class which implements the geometry of a tetrahedral element...
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
int64_t NNodes() const
Number of nodes of the mesh.
Definition: pzgmesh.h:126
virtual MElementType Type() const =0
Returns the element type acording to pzeltype.h.
TPZAdmChunkVector< TPZGeoNode > & NodeVec()
Definition: pzgmesh.h:140
Contains declaration of TPZGeoElRefPattern class which implements a generic geometric element which i...
TPZGeoMesh * fComputedGeomesh
A geometric mesh being computed.
Contains the TPZRefPattern class which defines the topology of the current refinement pattern to a me...
TPZGeoNode & Node(int i) const
Returns the ith node of the element.
Definition: pzgeoel.cpp:2570
Contains the TPZGeoCube class which implements the geometry of hexahedra element. ...
Contains the TPZGeoPoint class which implements the geometry of a point element or 0-D element...
virtual int NNodes() const =0
Returns the number of nodes of the element.
void BuildConnectivity()
Build the connectivity of the grid.
Definition: pzgmesh.cpp:967
TPZVec< TPZAutoPointer< TPZGeoMesh > > fSubBases
Vector of n bases to be connected to the original base.
int Dimension()
Get Dimension.
Definition: pzgmesh.h:291
bool fIsPrism
3d extrusion is based on prisms
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
int Id() const
Returns the identity of the current node.
Definition: pzgnode.h:68
bool fIsTetrahedron
3d extrusion is based on tetrahedrons
int ffrontMatID
Extrusion front material id.
Contains the TPZGeoPrism class which implements the geometry of a prism element.
TPZHierarquicalGrid & operator=(const TPZHierarquicalGrid &other)
bool fNonAffineQ
Non affine 2D and 3D extrusion.
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138