NeoPZ
tpzquadraticprism.cpp
Go to the documentation of this file.
1 
5 #include "tpzquadraticprism.h"
6 #include "tpzgeoblend.h"
7 #include "tpzgeoelmapped.h"
8 
9 #include "pzgeoelrefless.h.h"
10 #include "tpzgeoelrefpattern.h.h"
11 #include "pznoderep.h.h"
12 #include "pzshapepiram.h"
13 #include "tpzgeomid.h"
14 
15 #include "pzlog.h"
16 
17 #ifdef LOG4CXX
18 static LoggerPtr logger(Logger::getLogger("pz.specialmaps.quadraticprism"));
19 #endif
20 
21 #ifdef _AUTODIFF
22 #include "fad.h"
23 #endif
24 
25 using namespace pzshape;
26 using namespace pzgeom;
27 using namespace pztopology;
28 
29 template<class T>
30 void TPZQuadraticPrism::TShape(const TPZVec<T> &par,TPZFMatrix<T> &phi,TPZFMatrix<T> &dphi) {
31 
32  T qsi = par[0], eta = par[1], zeta = par[2];
33 
34  phi(0,0) = -0.5*(-1. + eta + qsi)*(-1. + zeta)*(2.*(eta + qsi) + zeta);
35  phi(1,0) = 0.5*qsi*(-1. + zeta)*(2. - 2.*qsi + zeta);
36  phi(2,0) = 0.5*eta*(-1. + zeta)*(2. - 2.*eta + zeta);
37  phi(3,0) = -0.5*(-1. + eta + qsi)*(1. + zeta)*(-2.*(eta + qsi) + zeta);
38  phi(4,0) = 0.5*qsi*(1. + zeta)*(-2. + 2.*qsi + zeta);
39  phi(5,0) = 0.5*eta*(1. + zeta)*(-2. + 2.*eta + zeta);
40  phi(6,0) = 2.*qsi*(-1. + eta + qsi)*(-1. + zeta);
41  phi(7,0) = -2.*eta*qsi*(-1. + zeta);
42  phi(8,0) = 2.*eta*(-1. + eta + qsi)*(-1. + zeta);
43  phi(9,0) = (-1. + eta + qsi)*(-1. + zeta*zeta);
44  phi(10,0) = qsi - qsi*zeta*zeta;
45  phi(11,0) = eta - eta*zeta*zeta;
46  phi(12,0) = -2.*qsi*(-1. + eta + qsi)*(1. + zeta);
47  phi(13,0) = 2.*eta*qsi*(1. + zeta);
48  phi(14,0) = -2.*eta*(-1. + eta + qsi)*(1. + zeta);
49 
50  dphi(0,0) = (1. - 2.*eta - 2.*qsi - 0.5*zeta)*(-1. + zeta);
51  dphi(1,0) = (1. - 2.*eta - 2.*qsi - 0.5*zeta)*(-1. + zeta);
52  dphi(2,0) = (-1. + eta + qsi)*(0.5 - eta - qsi - zeta);
53 
54  dphi(0,1) = -1. + qsi*(2. - 2.*zeta) + 0.5*zeta + 0.5*zeta*zeta;
55  dphi(1,1) = 0.;
56  dphi(2,1) = qsi*(0.5 - 1.*qsi + 1.*zeta);
57 
58  dphi(0,2) = 0.;
59  dphi(1,2) = -1. + eta*(2. - 2.*zeta) + 0.5*zeta + 0.5*zeta*zeta;
60  dphi(2,2) = eta*(0.5 - eta + zeta);
61 
62  dphi(0,3) = (-1. + 2.*eta + 2.*qsi - 0.5*zeta)*(1. + zeta);
63  dphi(1,3) = (-1. + 2.*eta + 2.*qsi - 0.5*zeta)*(1. + zeta);
64  dphi(2,3) = (-1. + eta + qsi)*(-0.5 + eta + qsi - zeta);
65 
66  dphi(0,4) = (-1. + 2.*qsi + 0.5*zeta)*(1. + zeta);
67  dphi(1,4) = 0.;
68  dphi(2,4) = qsi*(-0.5 + qsi + zeta);
69 
70  dphi(0,5) = 0.;
71  dphi(1,5) = -1. - 0.5*zeta + 0.5*zeta*zeta + eta*(2. + 2.*zeta);
72  dphi(2,5) = eta*(-0.5 + eta + zeta);
73 
74  dphi(0,6) = (-2. + 2.*eta + 4.*qsi)*(-1. + zeta);
75  dphi(1,6) = 2.*qsi*(-1. + zeta);
76  dphi(2,6) = 2.*qsi*(-1. + eta + qsi);
77 
78  dphi(0,7) = -2.*eta*(-1. + zeta);
79  dphi(1,7) = -2.*qsi*(-1. + zeta);
80  dphi(2,7) = -2.*eta*qsi;
81 
82  dphi(0,8) = 2.*eta*(-1. + zeta);
83  dphi(1,8) = (-2. + 4.*eta + 2.*qsi)*(-1. + zeta);
84  dphi(2,8) = 2.*eta*(-1. + eta + qsi);
85 
86  dphi(0,9) = -1. + zeta*zeta;
87  dphi(1,9) = -1. + zeta*zeta;
88  dphi(2,9) = 2.*(-1. + eta + qsi)*zeta;
89 
90  dphi(0,10) = 1. - zeta*zeta;
91  dphi(1,10) = 0.;
92  dphi(2,10) = -2.*qsi*zeta;
93 
94  dphi(0,11) = 0.;
95  dphi(1,11) = 1. - zeta*zeta;
96  dphi(2,11) = -2.*eta*zeta;
97 
98  dphi(0,12) = (2. - 2.*eta - 4.*qsi)*(1. + zeta);
99  dphi(1,12) = -2.*qsi*(1. + zeta);
100  dphi(2,12) = -2.*qsi*(-1. + eta + qsi);
101 
102  dphi(0,13) = 2.*eta*(1. + zeta);
103  dphi(1,13) = 2.*qsi*(1. + zeta);
104  dphi(2,13) = 2.*eta*qsi;
105 
106  dphi(0,14) = -2.*eta*(1. + zeta);
107  dphi(1,14) = (2. - 4.*eta - 2.*qsi)*(1. + zeta);
108  dphi(2,14) = -2.*eta*(-1. + eta + qsi);
109 }
110 
111 template<class T>
112 void TPZQuadraticPrism::X(const TPZFMatrix<REAL> &nodes,TPZVec<T> &loc,TPZVec<T> &x){
113 
114  TPZFNMatrix<15,T> phi(NNodes,1);
115  TPZFNMatrix<45,T> dphi(3,NNodes);
116  TShape(loc,phi,dphi);
117  int space = nodes.Rows();
118 
119  for(int i = 0; i < space; i++) {
120  x[i] = 0.0;
121  for(int j = 0; j < NNodes; j++) {
122  x[i] += phi(j,0)*nodes.GetVal(i,j);
123  }
124  }
125 
126 }
127 
128 template<class T>
129 void TPZQuadraticPrism::GradX(const TPZFMatrix<REAL> &nodes,TPZVec<T> &loc, TPZFMatrix<T> &gradx){
130 
131  gradx.Resize(3,3);
132  gradx.Zero();
133  int nrow = nodes.Rows();
134  int ncol = nodes.Cols();
135 #ifdef PZDEBUG
136  if(nrow != 3 || ncol != 15){
137  std::cout << "Objects of incompatible lengths, gradient cannot be computed." << std::endl;
138  std::cout << "nodes matrix must be 3x15." << std::endl;
139  DebugStop();
140  }
141 
142 #endif
143  TPZFNMatrix<15,T> phi(NNodes,1);
144  TPZFNMatrix<45,T> dphi(3,NNodes);
145  TShape(loc,phi,dphi);
146  for(int i = 0; i < NNodes; i++)
147  {
148  for(int j = 0; j < 3; j++)
149  {
150  gradx(j,0) += nodes.GetVal(j,i)*dphi(0,i);
151  gradx(j,1) += nodes.GetVal(j,i)*dphi(1,i);
152  gradx(j,2) += nodes.GetVal(j,i)*dphi(2,i);
153  }
154  }
155 
156 }
157 
162 // TPZGeoEl *TPZQuadraticPrism::CreateGeoElement(TPZGeoMesh &mesh, MElementType type,
163 // TPZVec<int64_t>& nodeindexes,
164 // int matid,
165 // int64_t& index)
166 // {
167 // return CreateGeoElementMapped(mesh,type,nodeindexes,matid,index);
168 // }
169 
170 // TPZGeoEl *TPZQuadraticPrism::CreateBCGeoEl(TPZGeoEl *orig,int side,int bc)
171 // {
172 // int ns = orig->NSideNodes(side);
173 // TPZManVector<int64_t> nodeindices(ns);
174 // int in;
175 // for(in=0; in<ns; in++)
176 // {
177 // nodeindices[in] = orig->SideNodeIndex(side,in);
178 // }
179 // int64_t index;
180 
181 // TPZGeoMesh *mesh = orig->Mesh();
182 // MElementType type = orig->Type(side);
183 
184 // TPZGeoEl *newel = mesh->CreateGeoBlendElement(type, nodeindices, bc, index);
185 // TPZGeoElSide me(orig,side);
186 // TPZGeoElSide newelside(newel,newel->NSides()-1);
187 
188 // newelside.InsertConnectivity(me);
189 // newel->Initialize();
190 
191 // return newel;
192 // }
193 
195 /* @param gmesh mesh in which the element should be inserted
196  @param matid material id of the element
197  @param lowercorner (in/out) on input lower corner o the cube where the element should be created, on exit position of the next cube
198  @param size (in) size of space where the element should be created
199  */
200 #include "tpzchangeel.h"
201 
202 void TPZQuadraticPrism::InsertExampleElement(TPZGeoMesh &gmesh, int matid, TPZVec<REAL> &lowercorner, TPZVec<REAL> &size)
203 {
204  TPZManVector<REAL,3> co(3),shift(3),scale(3);
205  TPZManVector<int64_t,4> nodeindexes(NCornerNodes);
206  for (int i=0; i<3; i++) {
207  scale[i] = size[i]/3.;
208  shift[i] = size[i]/2.+lowercorner[i];
209  }
210 
211  for (int i=0; i<NCornerNodes; i++) {
212  ParametricDomainNodeCoord(i, co);
213  co.Resize(3,0.);
214  for (int j=0; j<3; j++) {
215  co[j] = shift[j]+scale[j]*co[j]+(rand()*0.2/RAND_MAX)-0.1;
216  }
217  nodeindexes[i] = gmesh.NodeVec().AllocateNewElement();
218  gmesh.NodeVec()[nodeindexes[i]].Initialize(co, gmesh);
219  }
220  int64_t index;
221  gmesh.CreateGeoElement(EPrisma, nodeindexes, matid, index);
222  TPZGeoEl *gel = gmesh.Element(index);
223  int nsides = gel->NSides();
224  for (int is=0; is<nsides; is++) {
225  gel->SetSideDefined(is);
226  }
227  gel = TPZChangeEl::ChangeToQuadratic(&gmesh, index);
228  for (int node = gel->NCornerNodes(); node < gel->NNodes(); node++) {
229  TPZManVector<REAL,3> co(3);
230  gel->NodePtr(node)->GetCoordinates(co);
231  for (int i=0; i<3; i++) {
232  co[i] += (0.2*rand())/RAND_MAX - 0.1;
233  }
234  gel->NodePtr(node)->SetCoord(co);
235  }
236 }
237 
238 //void TPZQuadraticPrism::ParametricDomainNodeCoord(int node, TPZVec<REAL> &nodeCoord)
239 //{
240 // if(node > this->NNodes)
241 // {
242 // DebugStop();
243 // }
244 // nodeCoord.Resize(Dimension, 0.);
245 // switch (node) {
246 // case (0):
247 // {
248 // nodeCoord[0] = 0.;
249 // nodeCoord[1] = 0.;
250 // nodeCoord[2] = -1.;
251 // break;
252 // }
253 // case (1):
254 // {
255 // nodeCoord[0] = 1.;
256 // nodeCoord[1] = 0.;
257 // nodeCoord[2] = -1.;
258 // break;
259 // }
260 // case (2):
261 // {
262 // nodeCoord[0] = 0.;
263 // nodeCoord[1] = 1.;
264 // nodeCoord[2] = -1.;
265 // break;
266 // }
267 // case (3):
268 // {
269 // nodeCoord[0] = 0.;
270 // nodeCoord[1] = 0.;
271 // nodeCoord[2] = 1.;
272 // break;
273 // }
274 // case (4):
275 // {
276 // nodeCoord[0] = 1.;
277 // nodeCoord[1] = 0.;
278 // nodeCoord[2] = 1.;
279 // break;
280 // }
281 // case (5):
282 // {
283 // nodeCoord[0] = 0.;
284 // nodeCoord[1] = 1.;
285 // nodeCoord[2] = 1.;
286 // break;
287 // }
288 // case (6):
289 // {
290 // nodeCoord[0] = 0.5;
291 // nodeCoord[1] = 0.0;
292 // nodeCoord[2] = -1.0;
293 // break;
294 // }
295 // case (7):
296 // {
297 // nodeCoord[0] = 0.5;
298 // nodeCoord[1] = 0.5;
299 // nodeCoord[2] = -1.0;
300 // break;
301 // }
302 // case (8):
303 // {
304 // nodeCoord[0] = 0.0;
305 // nodeCoord[1] = 0.5;
306 // nodeCoord[2] = -1.0;
307 // break;
308 // }
309 // case (9):
310 // {
311 // nodeCoord[0] = 0.;
312 // nodeCoord[1] = 0.;
313 // nodeCoord[2] = 0.;
314 // break;
315 // }
316 // case (10):
317 // {
318 // nodeCoord[0] = 1.;
319 // nodeCoord[1] = 0.;
320 // nodeCoord[2] = 0.;
321 // break;
322 // }
323 // case (11):
324 // {
325 // nodeCoord[0] = 0.;
326 // nodeCoord[1] = 1.;
327 // nodeCoord[2] = 0.;
328 // break;
329 // }
330 // case (12):
331 // {
332 // nodeCoord[0] = 0.5;
333 // nodeCoord[1] = 0.0;
334 // nodeCoord[2] = 1.0;
335 // break;
336 // }
337 // case (13):
338 // {
339 // nodeCoord[0] = 0.5;
340 // nodeCoord[1] = 0.5;
341 // nodeCoord[2] = 1.0;
342 // break;
343 // }
344 // case (14):
345 // {
346 // nodeCoord[0] = 0.0;
347 // nodeCoord[1] = 0.5;
348 // nodeCoord[2] = 1.0;
349 // break;
350 // }
351 // default:
352 // {
353 // DebugStop();
354 // break;
355 // }
356 // }
357 //}
358 
360 
361 int TPZQuadraticPrism::ClassId() const{
362  return Hash("TPZQuadraticPrism") ^ TPZNodeRep<15,pztopology::TPZPrism>::ClassId() << 1;
363 }
364 
365 template class TPZRestoreClass< TPZGeoElRefPattern<TPZQuadraticPrism>>;
366 
367 /*@orlandini : I REALLY dont know why is this here, so I have commented the following lines.
368 If it breaks something, I am sorry.*/
369 //template class pzgeom::TPZNodeRep<15,TPZQuadraticPrism>;
370 
371 namespace pzgeom {
372  template void TPZQuadraticPrism::X(const TPZFMatrix<REAL>&, TPZVec<REAL>&, TPZVec<REAL>&);
373  template void TPZQuadraticPrism::GradX(const TPZFMatrix<REAL> &nodes,TPZVec<REAL> &loc, TPZFMatrix<REAL> &gradx);
374 
375 #ifdef _AUTODIFF
376  template void TPZQuadraticPrism::X(const TPZFMatrix<REAL>&, TPZVec<Fad<REAL> >&, TPZVec<Fad<REAL> >&);
377  template void TPZQuadraticPrism::GradX(const TPZFMatrix<REAL> &nodes,TPZVec<Fad<REAL> > &loc, TPZFMatrix<Fad<REAL> > &gradx);
378 #endif
379 
380 }
tpzquadraticprism.h
Contains the TPZQuadraticPrism class which defines a prism geometric element with quadratic map...
TPZAdmChunkVector::AllocateNewElement
int AllocateNewElement()
Makes more room for new elements.
Definition: pzadmchunk.h:184
TPZGeoNode::SetCoord
void SetCoord(const TPZVec< REAL > &x)
Sets all coordinates into the current node. It gets the dim values from x.
Definition: pzgnode.cpp:60
TPZGeoEl::NodePtr
TPZGeoNode * NodePtr(int i) const
Returns a pointer to the ith node of the element.
Definition: pzgeoel.cpp:2566
tpzchangeel.h
Contains the TPZChangeEl class. It is a special map.
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZGeoMesh::CreateGeoElement
virtual TPZGeoEl * CreateGeoElement(MElementType type, TPZVec< int64_t > &cornerindexes, int matid, int64_t &index, int reftype=1)
Generic method for creating a geometric element. Putting this method centrally facilitates the modifi...
Definition: pzgmesh.cpp:1296
pzshape
groups all classes dedicated to the computation of shape functions
Definition: pzshapeextend.h:16
TPZGeoEl::NCornerNodes
virtual int NCornerNodes() const =0
Returns the number of corner nodes of the element.
Fad
Definition: fad.h:54
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZGeoEl::NSides
virtual int NSides() const =0
Returns the number of connectivities of the element.
EPrisma
Definition: pzeltype.h:60
tpzgeoelmapped.h
Contains declaration of TPZGeoElMapped class which implements a geometric element using its ancestral...
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
pztopology
Groups all classes defining the structure of the master element.
Definition: PrismExtend.cpp:15
TPZChangeEl::ChangeToQuadratic
static TPZGeoEl * ChangeToQuadratic(TPZGeoMesh *Mesh, int64_t ElemIndex)
Turns an linear geoelement to quadratic.
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
tpzgeoblend.h
Contains the TPZGeoBlend class which implements a blending map from curved boundaries to the interior...
TPZGeoMesh::Element
TPZGeoEl * Element(int64_t iel)
Definition: pzgmesh.h:139
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZGeoMesh::NodeVec
TPZAdmChunkVector< TPZGeoNode > & NodeVec()
Definition: pzgmesh.h:140
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
co
REAL co[8][3]
Coordinates of the eight nodes.
Definition: tpzgensubstruct.cpp:56
pzgeom::TPZNodeRep
Implements ... Geometry Topology.
Definition: pznoderep.h:40
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZGeoEl::NNodes
virtual int NNodes() const =0
Returns the number of nodes of the element.
TPZGeoEl::SetSideDefined
virtual void SetSideDefined(int side)=0
Flags the side as defined, this means no neighbouring element was found.
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
pznoderep.h.h
Contains the implementation of the TPZNodeRep methods.
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
pzshapepiram.h
Contains TPZShapePiram class which implements the shape functions of a pyramid element.
tpzgeoelrefpattern.h.h
Contains the implementation of the TPZGeoElRefPattern methods.
TPZGeoNode::GetCoordinates
void GetCoordinates(TPZVec< REAL > &co)
Fill the coordinates of the node.
Definition: pzgnode.cpp:83
pzgeom
Groups all classes which model the geometry.
Definition: pzgeopoint.cpp:18
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
Generated on Sun Aug 16 2020 11:17:10 for NeoPZ by   doxygen 1.8.13