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