NeoPZ
pzgeoel.cpp
Go to the documentation of this file.
1 
6 #include "pzgeoel.h"
7 #include "pzcompel.h"
8 #include "pzgmesh.h"
9 #include "pzgnode.h"
10 #include "pzerror.h"
11 #include "pzmatrix.h"
12 #include "pzfmatrix.h"
13 #include "pztrnsform.h"
14 #include "pzvec.h"
15 #include "pzstack.h"
16 #include "pzquad.h"
17 #include "pzshapequad.h"
18 #include "pzshapetriang.h"
19 #include "TPZRefPattern.h"
20 #include "pzlog.h"
21 
22 #include "TPZGeoLinear.h"
23 #include "pzgeotriangle.h"
24 #include "pzgeoquad.h"
25 #include "pzgeotetrahedra.h"
26 #include "pzgeopyramid.h"
27 #include "pzgeoprism.h"
28 #include "TPZGeoCube.h"
29 #include "pzcmesh.h"
30 
31 
32 #include <stdio.h>
33 #include <stdlib.h>
34 
35 using namespace std;
36 using namespace pzgeom;
37 
38 #ifdef LOG4CXX
39 static LoggerPtr logger(Logger::getLogger("pz.mesh.tpzgeoel"));
40 static LoggerPtr loggerorient(Logger::getLogger("pz.mesh.tpzgeoel.orient"));
41 #endif
42 
43 
44 class TPZRefPattern;
45 
46 TPZFMatrix<REAL> TPZGeoEl::gGlobalAxes;
47 
48 // Destructor and Constructors
49 TPZGeoEl::~TPZGeoEl(){
50  int64_t index = Index();
51  if (this->fFatherIndex != -1) {
52  if(!this->Father()){
53  //Why did this element lose its father?
54  DebugStop();
55  } else {
56  int subelindex = WhichSubel();
57  if (subelindex == -1) {
58  DebugStop();
59  }
60  Father()->SetSubElement(subelindex, 0);
61  }
62  }
63  if (index != -1){
64  fMesh->ElementVec()[index] = NULL;
65  fMesh->ElementVec().SetFree(index); //the same line in TPZGeoMesh::DeleteElement was commented. Just call this once.
66  }
67  else
68  {
69  std::cout << __PRETTY_FUNCTION__ << " a derived class did reset the index of the element?\n";
70  }
71  fIndex = -1;
72  fId = -1;
73  fMatId = 0;
74  fNumInterfaces = 0;
75  fReference = 0;
76  fMesh = 0;
77 }
78 
79 
80 TPZGeoEl::TPZGeoEl(int64_t id,int materialid,TPZGeoMesh &mesh) {
81  fMesh = &mesh;
82  fId = id;
83  mesh.SetElementIdUsed(id);
84  fMatId = materialid;
85  this->fReference = NULL;
86  fFatherIndex = -1;
87  int64_t index = fMesh->ElementVec().AllocateNewElement();
88  fIndex = index;
89  fMesh->ElementVec()[index] = this;
90  this->fNumInterfaces = 0;
91 }
92 
93 TPZGeoEl::TPZGeoEl(const TPZGeoEl &el):TPZSavable(el){
94  this->fMesh = el.fMesh;
95  this->fId = fMesh->CreateUniqueElementId();
96  this->fMatId = el.fMatId;
97  this->fReference = NULL;
98  this->fFatherIndex = el.fFatherIndex;
99  this->fIndex = fMesh->ElementVec().AllocateNewElement();
100  this->fMesh->ElementVec()[fIndex] = this;
101  this->fNumInterfaces = 0;
102 }
103 
104 TPZGeoEl::TPZGeoEl(int materialid,TPZGeoMesh &mesh, int64_t &index) {
105  this->fMesh = &mesh;
106  this->fId = fMesh->CreateUniqueElementId();
107  this->fMatId = materialid;
108  this->fReference = NULL;
109  this->fFatherIndex = -1;
110  index = this->fMesh->ElementVec().AllocateNewElement();
111  this->fIndex = index;
112  this->fMesh->ElementVec()[fIndex] = this;
113  this->fNumInterfaces = 0;
114 }
115 
116 TPZGeoEl::TPZGeoEl(int materialid,TPZGeoMesh &mesh) {
117  this->fMesh = &mesh;
118  this->fId = fMesh->CreateUniqueElementId();
119  this->fMatId = materialid;
120  this->fReference = NULL;
121  this->fFatherIndex = -1;
122  const int64_t index = this->Mesh()->ElementVec().AllocateNewElement();
123  this->fIndex = index;
124  this->Mesh()->ElementVec()[index] = this;
125  this->fNumInterfaces = 0;
126 }
127 
128 void TPZGeoEl::Shape1d(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
129  if(num != 2 && num != 3){
130  PZError << "elcalc1d.shape, at this point only linear and quadratic elements\n";
131  return;
132  }
133 
134  if(num == 2) {
135  phi(0,0) = (1-x)/2.;
136  phi(1,0) = (1+x)/2.;
137  dphi(0,0) = -0.5;
138  dphi(0,1) = 0.5;
139  } else {
140  phi(0,0) = -x*(1.-x)*0.5;
141  dphi(0,0) = x-0.5;
142  phi(1,0) = (1.-x*x);
143  dphi(0,1) = -2.*x;
144  phi(2,0) = 0.5*x*(1.+x);
145  dphi(0,2) = x+0.5;
146  }
147 }
148 
149 void TPZGeoEl::ShapePhi1d(REAL x,int num,TPZFMatrix<REAL> &phi) {
150  if(num != 2 && num != 3){
151  PZError << "TPZGeoEl ShapePhi1d, at this point only linear and quadratic elements\n";
152  return;
153  }
154 
155  if(num == 2) {
156  phi(0,0) = (1-x)/2.;
157  phi(1,0) = (1+x)/2.;
158  } else {
159  phi(0,0) = -x*(1.-x)*0.5;
160  phi(1,0) = (1.-x*x);
161  phi(2,0) = 0.5*x*(1.+x);
162  }
163 }
164 
165 int TPZGeoEl::WhichSide(TPZVec<int64_t> &SideNodeIds) {
166  int64_t cap = SideNodeIds.NElements();
167  int nums = NSides();
168  for(int side=0; side<nums; side++) {
169  if(NSideNodes(side)==2 && cap == 2) {
170  int64_t isn1 = SideNodeIndex(side,0);
171  int64_t isn2 = SideNodeIndex(side,1);//sao = para side<3
172  if((isn1 == SideNodeIds[0] && isn2 == SideNodeIds[1]) ||
173  (isn2 == SideNodeIds[0] && isn1 == SideNodeIds[1])) return side;
174  } else if(NSideNodes(side)== 1 && cap ==1) {
175  if(SideNodeIndex(side,0) == SideNodeIds[0]) return side;
176  //completar
177  } else if(NSideNodes(side) == 3 && cap==3) {
178  int64_t sni[3],snx[3],k;
179  for(k=0;k<3;k++) snx[k] = SideNodeIndex(side,k);//el atual
180  for(k=0;k<3;k++) sni[k] = SideNodeIds[k];//el viz
181  for(k=0;k<3;k++) {
182  if(snx[0]==sni[k] && snx[1]==sni[(k+1)%3] && snx[2]==sni[(k+2)%3]) return side;
183  if(snx[0]==sni[k] && snx[1]==sni[(k+2)%3] && snx[2]==sni[(k+1)%3]) return side;
184  }//012 120 201 , 021 102 210
185  } else if(NSideNodes(side) == 4 && cap == 4) {//face quadrilateral
186  int64_t sni[4],snx[4],k;
187  for(k=0;k<4;k++) snx[k] = SideNodeIndex(side,k);//el atual
188  for(k=0;k<4;k++) sni[k] = SideNodeIds[k];//vizinho
189  if(snx[0]==sni[0]) {
190  for(k=1;k<4;k++) {
191  if(snx[1]==sni[k] && snx[2]==sni[k%3+1] && snx[3]==sni[(k+1)%3+1]) return side;
192  if(snx[1]==sni[k] && snx[2]==sni[(k+1)%3+1] && snx[3]==sni[k%3+1]) return side;
193  }// /* 0123 0231 0312 , 0132 0213 0321 */
194  } else if(snx[1]==sni[0]) {
195  for(k=1;k<4;k++) {
196  if(snx[0]==sni[k] && snx[2]==sni[k%3+1] && snx[3]==sni[(k+1)%3+1]) return side;
197  if(snx[0]==sni[k] && snx[2]==sni[(k+1)%3+1] && snx[3]==sni[k%3+1]) return side;
198  }// /* 0123 0231 0312 , 0132 0213 0321 */ /* 1023 1230 1302 , 1032 1203 1320 */
199  } else if(snx[2]==sni[0]) {
200  for(k=0;k<4;k++) {
201  if(snx[0]==sni[k] && snx[1]==sni[k%3+1] && snx[3]==sni[(k+1)%3+1]) return side;
202  if(snx[0]==sni[k] && snx[1]==sni[(k+1)%3+1] && snx[3]==sni[k%3+1]) return side;
203  }// /* 0123 0231 0312 , 0132 0213 0321 */ /* 2013 2130 2301 , 2031 2103 2310 */
204  } else if(snx[3]==sni[0]) {
205  for(k=0;k<4;k++) {
206  if(snx[0]==sni[k] && snx[1]==sni[k%3+1] && snx[2]==sni[(k+1)%3+1]) return side;
207  if(snx[0]==sni[k] && snx[1]==sni[(k+1)%3+1] && snx[2]==sni[k%3+1]) return side;
208  }// /* 0123 0231 0312 , 0132 0213 0321 */ / * 3012 3120 3201 , 3021 3102 3210 * /
209  }
210  } else if(cap<1 || cap > 4) {
211  int is;
212  for (is=0; is<nums; is++) {
213  if (NSideNodes(is) == cap) {
214  break;
215  }
216  }
217  if (is != nums) {
218  PZError << "TPZGeoEl::WhichSide must be extended\n";
219  DebugStop();
220  }
221  }
222  }
223  return -1;
224 }
225 
226 int TPZGeoEl::NeighbourExists(int side,const TPZGeoElSide &gel) {
227  TPZGeoElSide thisside(this,side);
228  if(gel == thisside) return 1;
229  TPZGeoElSide neighbour = Neighbour(side);
230  if(!neighbour.Exists()) return 0;
231  while(neighbour != thisside) {
232  if(gel == neighbour) return 1;
233  neighbour = neighbour.Neighbour();
234  }
235  return 0;
236 }
237 
238 void TPZGeoEl::Print(std::ostream & out) {
239 
240  out << "Element index " << fIndex << endl;
241  out << "Element id " << fId << endl;
242  out << "Is GeoBlend ? : ";
243  if(this->IsGeoBlendEl())
244  {
245  out << "true" << endl;
246  }
247  else
248  {
249  out << "false" << endl;
250  }
251  out << "Is Linear Mapping ? : ";
252  if(this->IsLinearMapping())
253  {
254  out << "true" << endl;
255  }
256  else
257  {
258  out << "false" << endl;
259  }
260  out << "Number of nodes " << NNodes() << endl;
261  out << "Corner nodes " << NCornerNodes() << endl;
262  out << "Nodes indexes ";
263  int i;
264  for (i = 0;i < NNodes();i++) out << NodeIndex(i) << " ";
265  out << "\nNumber of sides " << NSides() << endl;
266  if (fMatId < 0) out << "boundary condition " << fMatId << endl;
267  else out << "Material id " << fMatId << endl;
268  if (!Father()) out << "no father\n";
269  else out << "Father index " << Father()->Index() << endl;
270  if (!SubElement(0)) out << "no subelements";
271  else {
272  out << "Subelement indexes ";
273  for (i = 0;i < NSubElements();i++) {
274  if (!SubElement(i) ) continue;
275  out << SubElement(i)->Index() << ' ' ;
276  }
277  }
278  out << endl;
279  for (i = 0;i < NSides();i++) {
280  out << "Neighbours for side " << i << " : ";
281  TPZGeoElSide neighbour = Neighbour(i);
282  TPZGeoElSide thisside(this,i);
283  if (!neighbour.Exists())
284  {
285  out << "No neighbour\n";
286  }
287  else {
288  while (neighbour != thisside ) {
289  out << neighbour.Element()->Index() << "/" << neighbour.Side() << ' ';
290  neighbour = neighbour.Neighbour();
291  }
292  out << endl;
293  }
294  }
295  out << "Reference element pointer address: " << fReference << endl;
296  if (this->Reference()) out << "Reference element index : " << this->Reference()->Index() << endl;
297  out << "fNumInterfaces = " << fNumInterfaces << endl;
298 }
299 
300 void TPZGeoEl::PrintTopologicalInfo(std::ostream & out) {
301  int64_t elIndex = this->Index();
302  int nnodes = this->NNodes();
303  out << "Element index: " << elIndex << "\n";
304 
305  for(int n = 0; n < nnodes; n++) {
306  REAL nodeX = this->NodePtr(n)->Coord(0);
307  REAL nodeY = this->NodePtr(n)->Coord(1);
308  REAL nodeZ = this->NodePtr(n)->Coord(2);
309 
310  out << "Node " << n << " : " << nodeX << " , " << nodeY << " , " << nodeZ << "\n";
311  }
312 }
313 
314 std::ostream &operator<<(std::ostream &out,TPZGeoEl & el) {
315  el.Print(out);
316  return out;
317 }
318 
319 int TPZGeoEl::Level() {
320  int level = 0;
321  TPZGeoEl *father = Father();
322  while(father) {
323  level++;
324  father = father->Father();
325  }
326  return level;
327 }
328 
329 
330 int TPZGeoEl::GetTransformId2dQ(TPZVec<int> &idfrom,TPZVec<int> &idto) {
331 
332  if(idfrom[0]==idto[0] && idfrom[1]==idto[1]) return 0;//sentido horario : 0123
333  if(idfrom[0]==idto[0] && idfrom[1]==idto[3]) return 1;//sentido antihorario : 0321
334  //0 esta na posicao 3 e 1 esta na posicao 0 i.e. 1,2,3,0
335  if(idfrom[0]==idto[3] && idfrom[1]==idto[0]) return 2;//sentido horario : 1230
336  if(idfrom[0]==idto[1] && idfrom[1]==idto[0]) return 3;//sentido antihorario : 1032
337 
338  if(idfrom[0]==idto[2] && idfrom[1]==idto[3]) return 4;//sentido horario : 2301
339  if(idfrom[0]==idto[2] && idfrom[1]==idto[1]) return 5;//sentido antihorario : 2103
340 
341  if(idfrom[0]==idto[1] && idfrom[1]==idto[2]) return 6;//sentido horario : 3012
342  if(idfrom[0]==idto[3] && idfrom[1]==idto[2]) return 7;//sentido antihorario : 3210
343 
344  return 0;
345 }
346 
347 int TPZGeoEl::GetTransformId2dT(TPZVec<int> &idfrom,TPZVec<int> &idto) {
348  //REVISAR
349  if(idto[0]==idfrom[0] && idto[1]==idfrom[1]) return 0;//sentido horario
350  if(idto[0]==idfrom[0] && idto[1]==idfrom[2]) return 1;//sentido antihorario
351 
352  if(idto[0]==idfrom[1] && idto[1]==idfrom[2]) return 2;//sentido horario
353  if(idto[0]==idfrom[1] && idto[1]==idfrom[0]) return 3;//sentido antihorario
354 
355  if(idto[0]==idfrom[2] && idto[1]==idfrom[0]) return 4;//sentido horario
356  if(idto[0]==idfrom[2] && idto[1]==idfrom[1]) return 5;//sentido antihorario
357 
358  return 0;
359 }
360 
361 int TPZGeoEl::ElementExists(TPZGeoEl *elem,int64_t id) {
362 
363  if(elem == 0 && id == 0) {
364  PZError << "\nTPZGeoEl::ElementExists : element id or element pointer will not be null\n";
365  return -1;
366  }
367  int64_t nelg = fMesh->ElementVec().NElements(),index;
368  for(index=0;index<nelg;index++) {
369  TPZGeoEl *gel = fMesh->ElementVec()[index];
370  if(!gel) continue;
371  if(gel == elem || gel->Id() == id) return index;
372  }
373  return -1;
374 }
375 
376 TPZGeoElSide TPZGeoEl::Father2(int /*side*/) const{//Augusto:09/01/01
377  PZError << "TPZGeoEl::Father2 should never be called\n";
378  return TPZGeoElSide();
379 }
380 
381 int TPZGeoEl::FatherSide(int side, int son){
382  PZError << "TPZGeoEl::FatherSide should never be called\n";
383  return -1;
384 }
385 
386 TPZTransform<> TPZGeoEl::BuildTransform2(int /*side*/, TPZGeoEl * /*father*/, TPZTransform<> & /* tr */){//Augusto:09/01/01
387  PZError << "TPZGeoEl::BuildTransform2 should never be called\n";
388  return TPZTransform<>(0,0);
389 }
390 
391 
392 void TPZGeoEl::GetSubElements2(int /*side*/, TPZStack<TPZGeoElSide> &/*subel*/) const
393 {//Augusto:09/01/01
394  PZError << "TPZGeoEl::GetSubElements2 should never be called\n";
395 }
396 
397 void TPZGeoEl::GetSubElements2(int side, TPZStack<TPZGeoElSide> &subel, int dimension) const
398 {
399  TPZStack<TPZGeoElSide> subel2;
400  GetSubElements2(side,subel2);
401  int64_t cap = subel2.NElements();
402  int64_t s;
403  for(s=0; s<cap; s++) {
404  if(subel2[s].Dimension() == dimension) {
405  subel.Push(subel2[s]);
406  }
407  }
408 }
409 
410 void TPZGeoEl::GetAllSiblings(TPZStack<TPZGeoEl*> &unrefinedSons)
411 {
412 #ifdef PZDEBUG
413  PZError << __PRETTY_FUNCTION__ << "is deprecated. Use TPZGeoEl::YoungestChildren instead \n";
414  DebugStop();
415 #endif // PZDEBUG
416  int nsons = this->NSubElements();
417  for(int s = 0; s < nsons; s++)
418  {
419  TPZGeoEl * son = this->SubElement(s);
420  if(son->HasSubElement() == false)
421  {
422  unrefinedSons.Push(son);
423  }
424  else
425  {
426  son->GetAllSiblings(unrefinedSons);
427  }
428  }
429 }
430 
431 void TPZGeoEl::YoungestChildren(TPZStack<TPZGeoEl*> &unrefinedSons)
432 {
433  int nsons = this->NSubElements();
434  for(int s = 0; s < nsons; s++)
435  {
436  TPZGeoEl * son = this->SubElement(s);
437  if(son->HasSubElement() == false)
438  {
439  unrefinedSons.Push(son);
440  }
441  else
442  {
443  son->YoungestChildren(unrefinedSons);
444  }
445  }
446 }
447 
448 int TPZGeoEl::WhichSubel() const{
449 
450  if(fFatherIndex == -1) {
451  PZError << "TPZGeoEl::WhichSubel called with null element\n";
452  return -1;
453  }
454  int son;
455  TPZGeoEl *father = Father();
456  int nsub = father->NSubElements();
457  for(son=0;son<nsub;son++) if(father->SubElement(son) == this) break;
458  if(son > (nsub-1)){
459  PZError << "TPZGeoEl::WhichSubel son does not exist\n";
460 #ifdef LOG4CXX
461  {
462  std::stringstream sout;
463  sout << "Father element\n";
464  father->Print(sout);
465  sout << "Son element\n";
466  TPZGeoEl *foefel = (TPZGeoEl *) this;
467  foefel->Print(sout);
468  Mesh()->Print(sout);
469  LOGPZ_ERROR(logger,sout.str())
470  }
471 #endif
472  DebugStop();
473  return -1;
474  }
475  return son;
476 }
477 
478 int TPZGeoEl::WhichSide(TPZVec<REAL> &pt){
479  int dim = Dimension();
480  int nums = NSides();//n�mero de lados
481  REAL tol = 1.e-06;//toler�ncia do zero -> x � zero se -o<=x<=+o
482  int is;
483  for(is=0; is<nums; is++) {
484  int sdim = SideDimension(is);
485  TPZTransform<> t1 = SideToSideTransform(nums-1,is);
486  TPZTransform<> t2 = SideToSideTransform(is,nums-1);
487  TPZVec<REAL> pts(sdim),pt2(dim);
488  t1.Apply(pt,pts);
489  t2.Apply(pts,pt2);
490  REAL dif=0;
491  int d;
492  for(d=0; d<dim; d++) {
493  dif += (pt[d]-pt2[d])*(pt[d]-pt2[d]);
494  }
495  dif = sqrt(dif/dim);
496  if(dif < tol) return is;
497  }
498  cout << "TPZGeoEl::WhichSide ERROR : side not found" << endl ;
499 
500  return -1;
501 }
502 
503 void TPZGeoEl::CheckSubelDataStructure(){
504 
505  TPZVec<TPZGeoEl *> sub;
506  Divide(sub);
507  int side,nsubside,isub,ok; //no m�ximo 9 subelementos/lados
508  TPZStack<TPZGeoElSide> subside;//4 faces+4 arestas+1 canto
509  ofstream out("Checksubdata.dat");
510  for(side=0;side<NSides();side++){//percorre-se os lados do elemento atual = pai
511  GetSubElements2(side,subside);
512  nsubside = subside.NElements();
513  for(isub=0;isub<nsubside;isub++){
514  TPZGeoElSide neigh = subside[isub];//1o o subelemento depois os seus vizinhos
515  while(neigh.Element()){//existe vizinho!
516  TPZGeoElSide fatside = neigh.Element()->Father2(neigh.Side());
517  if(!fatside.Element()){
518  cout << "Inconsistencia de dados : \n";
519  cout << "Pai atual/lado : " << Id() << "/" << side << endl;
520  cout << "Pai do sub elemento e' nulo!\n";
521  cout << "Sub elemento/lado : " << neigh.Element()->Id() << "/" << neigh.Side() << endl;
522  cin >> ok;
523  } else //o pai existe � igual ao atual e a dimens�o dos lados � a mesma OK!
524  if(fatside.Element() == this){
525  cout << "Pai atual/lado : " << Id() << "/" << side << endl;
526  cout << "Sub elemento/lado : " << neigh.Element()->Id() << "/" << neigh.Side() << endl;
527  cout << "Pai sub elemento/lado : " << fatside.Element()->Id() << "/" << fatside.Side() << endl;
528  out << "Pai atual/lado : " << Id() << "/" << side << endl;
529  out << "Sub elemento/lado : " << neigh.Element()->Id() << "/" << neigh.Side() << endl;
530  out << "Pai sub elemento/lado : " << fatside.Element()->Id() << "/" << fatside.Side() << endl;
531  if(fatside.Side()!=side){
532  cout << "Dados acima inconsistentes : lados distintos\n";
533  out << "Dados acima inconsistentes : lados distintos\n";
534  cin >> ok;
535  }
536  } else {
537  cout << "OK!: Vizinho com pai diferente\n";
538  }
539  neigh = neigh.Neighbour();
540  if(!neigh.Element()){
541  cout << "Vizinho nulo\n\n";
542  out << "Vizinho nulo\n\n";
543  break;
544  }
545  cout << "\n";
546  out << "\n";
547  if(neigh.Element()==subside[isub].Element()) break;
548  //no pr�ximo while confia-se que o ciclo de conectividades � fechado e os viz sao filhos de alguem
549  TPZGeoEl *gel = neigh.Element()->Father2(neigh.Element()->NSides()-1).Element();//pai pelo interior:sempre existe
550  while(gel && !(gel== this)) neigh = neigh.Neighbour();
551  if(neigh.Element()==subside[isub].Element()){
552  cout << "Nao existe vizinho irmao!\n";
553  cout << "Sub elemento/lado " << neigh.Element()->Id() << "/" << neigh.Side() << endl;
554  cout << "OK!\n";
555  }
556  }
557  }
558  }
559  out.flush();
560  out.close();
561 }
562 
563 void TPZGeoEl::SetSubElementConnectivities() {
564 
565  int side;
566  for(side=0; side<NCornerNodes(); side++)
567  {
568  TPZGeoElSide thisside(this,side);
569  TPZStack<TPZGeoElSide> subel;
570  this->GetSubElements2(side,subel);
571  if(!subel[0].NeighbourExists(thisside))
572  {
573  subel[0].SetConnectivity(thisside);
574  }
575  }
576  for(side=NCornerNodes(); side<NSides(); side++)
577  {
578  TPZGeoElSide thisside(this,side);
579  TPZGeoElSide neighbour = this->Neighbour(side);
580  while(neighbour.Exists() && neighbour != thisside)
581  {
582  if(neighbour.HasSubElement() && neighbour.NSubElements() != 1) {
583  TPZStack<TPZGeoElSide> elvec,neighvec;
584  GetSubElements2(side,elvec);
585  neighbour.GetSubElements2(neighvec);
586 
587  // The currently divided element may have only one element as a son. In this case, the son is already
588  // neighbour of the father
589 
590  if(elvec.NElements() > 1) TPZGeoElSide::BuildConnectivities(elvec,neighvec);
591 
592  break;
593  }
594  neighbour = neighbour.Neighbour();
595  }
596  }
597  int nsubel = NSubElements();
598  for (int isub=0; isub<nsubel; isub++)
599  {
600  this->SubElement(isub)->InitializeNeighbours();
601  }
602 }
603 
604 REAL TPZGeoEl::CharacteristicSize()
605 {
606  int nn = NNodes();
607  if(!nn)
608  {
609  return 0.;
610  }
611  TPZVec<REAL> xmin(3),xmax(3);
612  TPZManVector<REAL,3> values(3);
613 
614  NodePtr(0)->GetCoordinates(values);
615  for(int c = 0; c < 3; c++)
616  {
617  xmin[c] = values[c];
618  xmax[c] = values[c];
619  }
620  for(int n = 1; n < nn; n++)
621  {
622  NodePtr(n)->GetCoordinates(values);
623  for(int c = 0; c < 3; c++)
624  {
625  xmin[c] = Min(values[c],xmin[c]);
626  xmax[c] = Max(values[c],xmax[c]);
627  }
628  }
629  REAL diagVecNorm = 0.;
630  for(int c = 0; c < 3; c++)
631  {
632  diagVecNorm += (xmax[c]-xmin[c])*(xmax[c]-xmin[c]);
633  }
634  diagVecNorm = sqrt(diagVecNorm);
635 
636  return diagVecNorm;
637 }
638 
639 REAL TPZGeoEl::SmallerEdge()
640 {
641  REAL norm = 0.;
642 
643  int firstEdge = this->NNodes();
644  TPZGeoElSide edge0(this,firstEdge);
645 
646  TPZVec<REAL> coords0(3), coords1(3);
647  int64_t node0id = edge0.SideNodeIndex(0);
648  int64_t node1id = edge0.SideNodeIndex(1);
649  this->Mesh()->NodeVec()[node0id].GetCoordinates(coords0);
650  this->Mesh()->NodeVec()[node1id].GetCoordinates(coords1);
651  for(int c = 0; c < 3; c++)
652  {
653  REAL delta = coords1[c] - coords0[c];
654  norm += delta*delta;
655  }
656  norm = sqrt(norm);
657 
658  for(int s = firstEdge+1; s < this->NSides(); s++)
659  {
660  TPZGeoElSide edgeOther(this,s);
661  if(edgeOther.Dimension() > 1)
662  {
663  break;
664  }
665 
667  node0id = edgeOther.SideNodeIndex(0);
668  node1id = edgeOther.SideNodeIndex(1);
669  this->Mesh()->NodeVec()[node0id].GetCoordinates(coords0);
670  this->Mesh()->NodeVec()[node1id].GetCoordinates(coords1);
671  REAL normTemp = 0.;
672  for(int c = 0; c < 3; c++)
673  {
674  REAL delta = coords1[c] - coords0[c];
675  normTemp += delta*delta;
676  }
677  normTemp = sqrt(normTemp);
678 
679  norm = Min(norm,normTemp);
681  }
682 
683  return norm;
684 }
685 
686 bool TPZGeoEl::ComputeXInverse(TPZVec<REAL> &XD, TPZVec<REAL> &qsi, REAL Tol) {
687  REAL error = 10.;
688  int iter = 0;
689  const int nMaxIter = 10000;
690  REAL radius = CharacteristicSize();
691  int dim = Dimension();
692  TPZManVector<REAL,3> X0(3);
693 // ZeroTolerance();
694  // First verify if the entry qsi yields the right point
695  if(qsi.NElements()!= dim)
696  {
697  PZError << "\nTPZGeoEl::ComputeXInverse vector dimension error\n";
698  qsi.Resize(Dimension(),0.);//zero esta em todos os elementos mestres
699  }
700  X(qsi,X0);//qsi deve ter dimensao do elemento atual
701  TPZFNMatrix<9> DelX(3,1);
702  int i;
703  for(i=0; i<3; i++)
704  {
705  DelX(i,0) = XD[i]-X0[i];
706  }
707  error = Norm(DelX)/radius;
708  if(error <= Tol)
709  {
710  return ( this->IsInParametricDomain(qsi) );
711  }
712 
713  // Verify if the point is not a corner node
714  TPZManVector<REAL,3> values(3);
715  int nn = NCornerNodes();
716  if(nn <= 1)
717  {
718  return false;
719  }
720  int in;
721  for(in=0; in<nn; in++)
722  {
723  NodePtr(in)->GetCoordinates(values);
724  TPZFNMatrix<3> DelX(3,1);
725  for(i=0; i<3; i++) DelX(i,0) = XD[i]-values[i];
726  error = Norm(DelX)/radius;
727  if(error <= Tol)
728  {
729  TPZVec<REAL> zero(0);
730  TPZTransform<> tr = SideToSideTransform(in, NSides()-1);
731  tr.Apply(zero, qsi);
732  return true;
733  }
734  }
735 
736  while(error > Tol && iter < nMaxIter)
737  {
738  iter++;
739  TPZFNMatrix<9> residual(dim,1),delqsi(dim,1);
740  REAL detJ;
741  TPZFNMatrix<9> J(dim,dim,0.),axes(dim,3,0.),Inv(dim,dim,0.);
742  TPZFNMatrix<9> JXt(dim,3,0.),JX(3,dim,0.),JXtJX(dim,dim,0.);
743  TPZFNMatrix<9,REAL> gradx(3,dim);
744  GradX(qsi, gradx);
745  Jacobian(gradx,J,axes,detJ,Inv);
746  if(fabs(detJ) < 2.e-10)
747  {
748  TPZManVector<REAL,3> center(Dimension(),0.);
749  CenterPoint(NSides()-1, center);
750  cout << "ComputeXInverse found zero Jacobian Index " << this->fIndex << " qsi " << qsi << " detJ " << detJ << std::endl;
751  for(int ik = 0; ik < qsi.NElements(); ik++)
752  {
753  qsi[ik] += Tol*(center[ik]-qsi[ik]);
754  }
755  residual(0,0) = 1.e12;
756  }
757  else
758  {
759  TPZFNMatrix<9> axest;
760  axes.Transpose(&axest);
761  axest.Resize(3,dim);//casos 1D e 2D onde JX espacial � 1x3 e 2x3 respectivamente
762  if(dim==1)
763  {
764  JX(0,0) = axest(0,0)*J(0,0);
765  JX(1,0) = axest(1,0)*J(0,0);
766  JX(2,0) = axest(2,0)*J(0,0);
767  }
768  else
769  {
770  axest.Multiply(J,JX,0);
771  }
772 
773  JX.Transpose(&JXt);
774  JXt.Multiply(JX,JXtJX,0);//JXtJX = JXt*JX;
775  JXt.Multiply(DelX,residual);//cout << "\nComputeXInverse: : \n";
776  JXtJX.SolveDirect(residual,ELU);//cout << "Atual/dimensao : " << Id() << " / " << Dimension();
777  for(i=0; i<dim; i++)
778  {
779  qsi[i] += residual(i,0);
780  }
781 
782  }
783  X(qsi,X0);
784  for(i=0; i<3; i++)
785  {
786  DelX(i,0) = XD[i]-X0[i];
787  }
788  //A norma sobre coordenada parametrica eh mais objetiva pois os limites do
789  //elemento mestre sao mais claros que os do elemento real
790  // error = Norm(DelX);
791  error = Norm(DelX)/radius;
792  }
793 
794 #ifdef PZDEBUG
795  if(iter == nMaxIter)
796  {
797  std::stringstream sout;
798  sout << "Error at " << __PRETTY_FUNCTION__ << " - nMaxIter was reached before tolerance is achieved - ElementId" << this->Id() << std::endl;
799  PZError << "\n" << sout.str() << "\n";
800  Print(std::cout);
801  int nnodes = NNodes();
802  for (int i=0; i<NNodes(); i++) {
803  NodePtr(i)->Print();
804  }
805 
806 #ifdef LOG4CXX
807  LOGPZ_ERROR(logger,sout.str().c_str());
808 #endif
809  }
810 #endif
811 
812  return ( this->IsInParametricDomain(qsi) );
813 }
814 
815 
816 void TPZGeoEl::TransformSonToFather(TPZGeoEl *ancestor, TPZVec<REAL> &qsiSon, TPZVec<REAL> &qsiAncestor){
817 
818  int dim = Dimension();
819  if(qsiAncestor.NElements()!= dim)
820  {
821  PZError << "\nTPZGeoEl::TransformSonToFather vector dimension error\n";
822  qsiAncestor.Resize(Dimension(),0.);//zero esta em todos os elementos mestres
823  }
824 
825  TPZVec<REAL> xson;
826  X(qsiSon,xson);
827 
828  TPZGeoEl *father = this;
829  while(father != ancestor)
830  {
831  father = father->Father();
832  }
833  if(!father)
834  {
835  DebugStop();
836  }
837  TPZTransform<> tr(dim);
838  tr = BuildTransform2(NSides()-1, father, tr);
839  tr.Apply(qsiSon, qsiAncestor);
840  REAL Tol;
841  ZeroTolerance(Tol);
842  father->ComputeXInverse(xson, qsiAncestor,Tol);
843 }
844 
845 TPZTransform<> TPZGeoEl::ComputeParamTrans(TPZGeoEl *fat,int fatside, int sideson){
846 
847  //transformacao do lado de elemento pequeno para elemento grande que o contem
848  int dimf = fat->Dimension();
849  int dimsf = fat->SideDimension(fatside);
850  int dim = Dimension();
851  int dimss = SideDimension(sideson);
852  if(dimsf < dimss){
853  PZError << "\nTPZGeoEl::ComputeParamTrans called with sides error\n";
854  DebugStop();
855  }
856 
858  if(!fat->SideDimension(fatside)) return TPZTransform<>(0,0);
859 
860  REAL weight;
861  TPZFNMatrix<9> jac(dim,dim),axes(3,3,0.);
862  TPZFNMatrix<9> jacinv(dim,dim);
863  TPZManVector<REAL,3> x(3,0.);
864  TPZManVector<REAL,3> intpoint(dimss,0.);
865  int tam = (dimss+1);
866  TPZFNMatrix<16> hess(tam,tam,0.),grad0(tam,1,0.);
867  TPZIntPoints *intrule = CreateSideIntegrationRule(sideson,2);
868  TPZManVector<int,3> order(dimss,2);
869  intrule->SetOrder(order);
870  //integra��o sobre o lado-filho contido no lado-pai
871  int ij,ik,indp;
872  REAL D2Edaikdaij,D2Edcidaij,D2Edci2;
873  D2Edcidaij = 0.;
874  D2Edci2 = 0.;
875  for(ij=0;ij<dimss;ij++){
876  for(ik=ij;ik<dimss;ik++){
877  D2Edaikdaij = 0.;
878  D2Edcidaij = 0.;
879  for(indp = 0; indp < intrule->NPoints(); ++indp){
880  intrule->Point(indp,intpoint,weight);
881  //Jacobian(intpoint,jac,axes,detjac,jacinv);/**estes passos n�o s�o precisos pois:*/
882  //weight *= fabs(detjac);/**ambas as matrizes, A e b, multiplicam a mesma constante detjac*/
883  D2Edaikdaij += intpoint[ik]*intpoint[ij]*weight;
884  if(ik==ij) D2Edcidaij += intpoint[ij]*weight;
885  if(ij==0 && ik==0) D2Edci2 += weight;
886  }
887  hess(ij,ik) = 2.*D2Edaikdaij;
888  hess(ik,ij) = hess(ij,ik);
889  if(ik==ij) {
890  hess(ij,dimss) = 2.*D2Edcidaij;
891  hess(dimss,ij) = hess(ij,dimss);
892  }
893  if(ij==0 && ik==0) hess(dimss,dimss) = 2.*D2Edci2;
894  }
895  }// final do integral hess
896  //do lado sideson para o elemento atual (filho)
897  TPZTransform<> tsidetoson(Dimension());//identidade
898  if(dimss<Dimension()) tsidetoson = SideToSideTransform(sideson,NSides()-1);
899  TPZTransform<> fatelside = fat->SideToSideTransform(fat->NSides()-1,fatside);
900  TPZManVector<REAL,3> sidepoint(Dimension());//dimensao do dominio da transformacao X do filho
901  int j;//transf. para o lado do pai
902  TPZFNMatrix<9> A(dimsf,dimss,0.),sol(dimsf,1,0.);
903  REAL Tol;
904  ZeroTolerance(Tol);
905  for(int ifat=0;ifat<dimsf;ifat++){//numero de variaveis do pai
906  REAL DEdci = 0.;
907  for(j=0;j<(dimss+1);j++){
908  REAL DEdaij = 0.;
909  for(indp = 0; indp < intrule->NPoints(); ++indp){
910  intrule->Point(indp,intpoint,weight);
911  tsidetoson.Apply(intpoint,sidepoint);//lado do filho para o seu interior: mestre
912  X(sidepoint,x);//ponto do mestre do filho para o filho deformado, 3 coordenadas
913  TPZVec<REAL> csi(dimf,0.);
914  fat->ComputeXInverse(x,csi,Tol);//csi ponto no pai mestre
915  TPZVec<REAL> outcsi(dimsf,0.);
916  fatelside.Apply(csi,outcsi);
917  // fat->ProjectPointElToSide(fatside,csi,outcsi);
918  if(j<dimss) DEdaij += outcsi[ifat]*intpoint[j]*weight;
919  if(j==0) DEdci += outcsi[ifat]*weight;
920  }
921  if(j<dimss) grad0(j,0) = 2.*DEdaij;
922  if(j==0) grad0(dimss,0) = 2.*DEdci;
923  if(!dimss) grad0(0,0) = DEdci;
924  }//final integral gradiente
925  //resolu��o do sistema para cada variavel ifat do pai
926  if(dimss) hess.SolveDirect(grad0,ELU);
927  for(int k=0;k<dimss;k++) A(ifat,k) = grad0(k,0);
928  sol(ifat,0) = grad0(dimss,0);
929  }//fim sistema ifat
930  delete intrule;
931  TPZTransform<> t(dimsf,dimss);
932  t.SetMatrix(A,sol);
933  return t;
934 }
935 
936 REAL TPZGeoEl::Distance(TPZVec<REAL> &centel,TPZVec<REAL> &centface){
937 
938  if(centel.NElements() != 3 || centface.NElements() != 3)
939  PZError << "TPZGeoEl::Distance point dimension error\n";
940 
941  int i;
942  REAL distance = 0.;
943  for(i=0;i<3;i++)
944  distance += (centel[i]-centface[i])*(centel[i]-centface[i]);
945 
946  return sqrt(distance);
947 }
948 
949 REAL TPZGeoEl::ElementRadius(){
950 
951  switch( this->Dimension() )
952  {
953  case 0:
954  return 0.;
955  // break;
956 
957  case 1:
958  case 2:{
959  TPZManVector<REAL, 3> centel (this->Dimension(), 0.);
960  TPZManVector<REAL, 3> centface(this->Dimension(), 0.);
961  TPZManVector<REAL, 3> masscent(3,0.), massface(3, 0.);
962  REAL mindist = 1000.;
963  REAL dist;
964 
965  CenterPoint(NSides()-1,centel);
966  X(centel,masscent);
967 
968  int nsides = NSides();
969  for(int iside = 0; iside < nsides - 1; iside++){
970  CenterPoint(iside, centface);
971  X(centface,massface);
972  dist = Distance(masscent,massface);
973  if(mindist > dist) mindist = dist;
974  }
975  return mindist;
976  // break;
977  }
978 
979  case 3:{
980  TPZManVector<REAL, 3> centel(3,0.),centface(3,0.);
981  TPZManVector<REAL, 3> masscent(3,0.),massface(3,0.);
982  REAL mindist = 1000.,dist;
983  CenterPoint(NSides()-1,centel);
984  X(centel,masscent);
985  int nfaces,face,face0;
986  if(NSides()==15) {nfaces = 14; face0 = 10;}//tetrahedron
987  else if(NSides()==19) {nfaces = 18; face0 = 13;}//pyramid
988  else if(NSides()==21) {nfaces = 20; face0 = 15;}//prism
989  else if(NSides()==27) {nfaces = 26; face0 = 20;}//hexahedro
990  else return 0.;//line, point, triangle, quadrilateral
991  for(face=face0;face<nfaces;face++){
992  CenterPoint(face,centface);
993  X(centface,massface);
994  dist = Distance(masscent,massface);
995  if(mindist > dist) mindist = dist;
996  }
997  return mindist;
998  // break;
999  }
1000 
1001  default:
1002  PZError << "TPZGeoEl::ElementRadius - Dimension not implemented." << endl;
1003  return 0.;
1004 
1005  }//end of switch
1006  // return 0.;
1007 } //end of method
1008 
1009 REAL TPZGeoEl::TriangleArea(TPZVec<TPZGeoNode *> &nodes){
1010 
1011  if(nodes.NElements() != 3 && nodes.NElements() != 4){
1012  PZError << "TPZGeoEl::AreaFromTheFaceT argument error size: nodes\n";
1013  nodes.Resize(0);
1014  return 0.;
1015  }
1016 
1017  REAL cb0 = nodes[2]->Coord(0) - nodes[1]->Coord(0);
1018  REAL cb1 = nodes[2]->Coord(1) - nodes[1]->Coord(1);
1019  REAL cb2 = nodes[2]->Coord(2) - nodes[1]->Coord(2);
1020  REAL ab0 = nodes[0]->Coord(0) - nodes[1]->Coord(0);
1021  REAL ab1 = nodes[0]->Coord(1) - nodes[1]->Coord(1);
1022  REAL ab2 = nodes[0]->Coord(2) - nodes[1]->Coord(2);
1023 
1024  //produto vetorial
1025  REAL coord0 = cb1*ab2-ab1*cb2;
1026  REAL coord1 = ab0*cb2-cb0*ab2;
1027  REAL coord2 = cb0*ab1-ab0*cb1;
1028  //norma da metade do vetor
1029  return ( 0.5*sqrt(coord0*coord0+coord1*coord1+coord2*coord2) );
1030 }
1031 
1032 REAL TPZGeoEl::QuadArea(TPZVec<TPZGeoNode *> &nodes){
1033  if(nodes.NElements() != 4){
1034  PZError << "TPZGeoEl::AreaFromTheFaceQ argument error size: nodes\n";
1035  nodes.Resize(0);
1036  return 0.;
1037  }
1038 
1039  REAL areat1 = TriangleArea(nodes);
1040 
1041  nodes[1] = nodes[0];
1042 
1043  nodes[0] = nodes[2];
1044  nodes[2] = nodes[1];//antigo node0
1045  nodes[1] = nodes[3];
1046 
1047  REAL areat2 = TriangleArea(nodes);
1048 
1049  return (areat1+areat2);
1050 }
1051 
1052 REAL TPZGeoEl::Volume(){
1053 
1054  TPZManVector<REAL,3> param(3,0.);
1055  REAL detjac;
1056  TPZFNMatrix<9> jacinv(3,3),jacobian(3,3),axes(3,3), gradx(3,3);
1057  //supondo jacobiano constante: X linear
1058  CenterPoint(NSides()-1,param);
1059  GradX(param, gradx);
1060  Jacobian(gradx,jacobian,axes,detjac,jacinv);
1061  return (RefElVolume()*detjac);//RefElVolume(): volume do elemento de refer�ncia
1062 }
1063 
1064 REAL TPZGeoEl::SideArea(int side){
1065 
1066  if(side < 0 || side > NSides()-1)
1067  PZError << "TPZGeoEl::AreaFromTheFace side error, side = " << side << endl;
1068 
1069  if(SideDimension(side) != 2)
1070  {
1071  TPZGeoElSide gelside(this,side);
1072  return gelside.Area();
1073  }
1074 
1075  if(SideDimension(side) == 2){
1076 
1077  int nsn = NSideNodes(side);
1078 
1079  TPZVec<TPZGeoNode *> nodes(nsn);
1080  int i;
1081 
1082  for(i=0;i</*3*/nsn;i++)
1083  nodes[i] = &Mesh()->NodeVec()[ SideNodeIndex(side,i) ];
1084  if (nsn==4)
1085  return ( QuadArea(nodes) );
1086  else
1087  return ( TriangleArea(nodes) );
1088  }
1089  return 0.;
1090 }
1091 
1092 TPZCompEl *TPZGeoEl::CreateBCCompEl(int side,int bc,TPZCompMesh &cmesh) {
1093  TPZGeoEl *gel = CreateBCGeoEl(side,bc);
1094  int64_t index;
1095  return cmesh.CreateCompEl(gel,index);
1096 }
1097 
1098 void TPZGeoEl::RemoveConnectivities(){
1099 
1100  int nsides = NSides(),side;
1101  for(side=0;side<nsides;side++){
1102  TPZGeoElSide thisside(this,side);
1103  TPZGeoElSide neighbour (thisside.Neighbour());
1104  thisside.RemoveConnectivity();
1105  if(neighbour != thisside){
1106  TPZGeoElSide neighneigh = neighbour;
1107  do{
1108  if(neighneigh.ResetBlendConnectivity(fIndex)){
1109  neighneigh.Element()->SetNeighbourForBlending(neighneigh.Side());
1110  }
1111  neighneigh = neighneigh.Neighbour();
1112  } while (neighbour != neighneigh);
1113  }
1114  }
1115 }
1116 
1117 void TPZGeoEl::InitializeNeighbours(){
1118  int i,j;
1119  for (i=0;i<NSides();i++){
1120 
1121  TPZStack <TPZGeoElSide> subel;
1122  if (HasSubElement()){
1123  GetSubElements2(i,subel);
1124  for (j=0;j<subel.NElements();j++){
1125  TPZGeoEl *el = subel[j].Element();
1126  el->InitializeNeighbours();
1127  }
1128  }
1129  TPZGeoElSide neighside = Neighbour(i);
1130  if (!neighside.Element() || neighside.Side() == -1){
1131  TPZGeoElSide thisside (this,i);
1132  SetNeighbour(i,thisside);
1133  }
1134  }
1135 }
1136 
1137 void TPZGeoEl::MidSideNodeIndices(int side,TPZVec<int64_t> &indices) const {
1138  indices.Resize(1);
1139  MidSideNodeIndex(side,indices[0]);
1140  if(indices[0] == -1) indices.Resize(0);
1141 }
1142 
1144 void TPZGeoEl::Jacobian(TPZVec<REAL> &qsi,TPZFMatrix<REAL> &jac,TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv) const{
1145  TPZFNMatrix<9,REAL> gradx;
1146  GradX(qsi, gradx);
1147 #ifdef LOG4CXX
1148  if(logger->isDebugEnabled())
1149  {
1150  std::stringstream sout;
1151  gradx.Print("gradx",sout);
1152  LOGPZ_DEBUG(logger, sout.str())
1153  }
1154 #endif
1155  Jacobian(gradx, jac, axes, detjac, jacinv);
1156 }
1157 
1159 void TPZGeoEl::JacobianXYZ(TPZVec<REAL> &qsi,TPZFMatrix<REAL> &jac,TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv) const{
1160  TPZFNMatrix<9,REAL> gradx;
1161  GradX(qsi, gradx);
1162  JacobianXYZ(gradx, jac, axes, detjac, jacinv);
1163 }
1164 
1165 
1166 void TPZGeoEl::Jacobian(const TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &jac,TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv){
1167 
1168  detjac = 0.0;
1169  int nrows = gradx.Rows();
1170  int ncols = gradx.Cols();
1171  int dim = ncols;
1172 
1173  switch (dim) {
1174  case 0:
1175  {
1176  jac.Resize(dim,dim);
1177  axes.Resize(dim,3);
1178  jacinv.Resize(dim,dim);
1179  detjac=1.;
1180 
1181  break;
1182  }
1183 
1184  case 1:
1185  {
1186  jac.Resize(dim,dim);
1187  axes.Resize(dim,3);
1188  jacinv.Resize(dim,dim);
1189  jac.Zero();
1190 
1192  TPZManVector<REAL,3> v_1(3,0.);
1193 
1194  for (int i = 0; i < nrows; i++) {
1195  v_1[i] = gradx.GetVal(i,0);
1196  }
1197 
1198  REAL norm_v_1 = 0.;
1199  for(int i = 0; i < nrows; i++) {
1200  norm_v_1 += v_1[i]*v_1[i];
1201  }
1202 
1203  norm_v_1 = sqrt(norm_v_1);
1204  jac(0,0) = norm_v_1;
1205  detjac = norm_v_1;
1206 
1207  if(IsZero(detjac))
1208  {
1209 
1210 #ifdef PZDEBUG
1211  std::stringstream sout;
1212  sout << "Singular Jacobian, 1 determinant of jacobian = " << detjac << std::endl;
1213  LOGPZ_ERROR(logger, sout.str())
1214  //DebugStop();
1215 #endif
1216  detjac = ZeroTolerance();
1217  }
1218 
1219  jacinv(0,0) = 1.0/detjac;
1220 
1221  for(int i=0; i < 3; i++) {
1222  axes(0,i) = v_1[i]/norm_v_1;
1223  }
1224 
1225  }
1226  break;
1227  case 2:
1228  {
1229 
1230  jac.Resize(dim,dim);
1231  axes.Resize(dim,3);
1232  jacinv.Resize(dim,dim);
1233  jac.Zero();
1234 
1236  TPZManVector<REAL,3> v_1(3,0.), v_2(3,0.);
1237 
1239  TPZManVector<REAL,3> v_1_til(3,0.), v_2_til(3,0.);
1240 
1241  for (int i = 0; i < nrows; i++) {
1242  v_1[i] = gradx.GetVal(i,0);
1243  v_2[i] = gradx.GetVal(i,1);
1244  }
1245 
1246  REAL norm_v_1_til = 0.0;
1247  REAL norm_v_2_til = 0.0;
1248  REAL v_1_dot_v_2 = 0.0;
1249 
1250  for(int i = 0; i < 3; i++) {
1251  norm_v_1_til += v_1[i]*v_1[i];
1252  v_1_dot_v_2 += v_1[i]*v_2[i];
1253  }
1254  norm_v_1_til = sqrt(norm_v_1_til);
1255 
1256  for(int i=0 ; i < 3; i++) {
1257  v_1_til[i] = v_1[i] / norm_v_1_til; // Normalizing
1258  v_2_til[i] = v_2[i] - v_1_dot_v_2 * v_1_til[i] / norm_v_1_til;
1259  norm_v_2_til += v_2_til[i]*v_2_til[i];
1260  }
1261  norm_v_2_til = sqrt(norm_v_2_til);
1262 
1263 
1264  jac(0,0) = norm_v_1_til;
1265  jac(0,1) = v_1_dot_v_2/norm_v_1_til;
1266  jac(1,1) = norm_v_2_til;
1267 
1268  detjac = jac(0,0)*jac(1,1)-jac(1,0)*jac(0,1);
1269 
1270  jacinv(0,0) = +jac(1,1)/detjac;
1271  jacinv(1,1) = +jac(0,0)/detjac;
1272  jacinv(0,1) = -jac(0,1)/detjac;
1273  jacinv(1,0) = -jac(1,0)/detjac;
1274 
1275  if(IsZero(detjac))
1276  {
1277 
1278 #ifdef PZDEBUG
1279  std::stringstream sout;
1280  sout << "Singular Jacobian, 2 determinant of jacobian = " << detjac << std::endl;
1281  LOGPZ_ERROR(logger, sout.str())
1282  DebugStop();
1283 #endif
1284  detjac = ZeroTolerance();
1285  }
1286 
1287  for(int i=0; i < 3; i++) {
1288  v_2_til[i] /= norm_v_2_til; // Normalizing
1289  axes(0,i) = v_1_til[i];
1290  axes(1,i) = v_2_til[i];
1291  }
1292 
1293  }
1294  break;
1295  case 3:
1296  {
1297  jac.Resize(dim,dim);
1298  axes.Resize(dim,3);
1299  jacinv.Resize(dim,dim);
1300  jac.Zero();
1301 
1302  for (int i = 0; i < nrows; i++) {
1303  jac(i,0) = gradx.GetVal(i,0);
1304  jac(i,1) = gradx.GetVal(i,1);
1305  jac(i,2) = gradx.GetVal(i,2);
1306  }
1307 
1308  detjac -= jac(0,2)*jac(1,1)*jac(2,0);//- a02 a11 a20
1309  detjac += jac(0,1)*jac(1,2)*jac(2,0);//+ a01 a12 a20
1310  detjac += jac(0,2)*jac(1,0)*jac(2,1);//+ a02 a10 a21
1311  detjac -= jac(0,0)*jac(1,2)*jac(2,1);//- a00 a12 a21
1312  detjac -= jac(0,1)*jac(1,0)*jac(2,2);//- a01 a10 a22
1313  detjac += jac(0,0)*jac(1,1)*jac(2,2);//+ a00 a11 a22
1314 
1315  if(IsZero(detjac))
1316  {
1317 
1318 #ifdef PZDEBUG
1319  std::stringstream sout;
1320  sout << "Singular Jacobian, 3 determinant of jacobian = " << detjac << std::endl;
1321  LOGPZ_ERROR(logger, sout.str())
1322  DebugStop();
1323 #endif
1324  detjac = ZeroTolerance();
1325  }
1326 
1327  jacinv(0,0) = (-jac(1,2)*jac(2,1)+jac(1,1)*jac(2,2))/detjac;//-a12 a21 + a11 a22
1328  jacinv(0,1) = ( jac(0,2)*jac(2,1)-jac(0,1)*jac(2,2))/detjac;//a02 a21 - a01 a22
1329  jacinv(0,2) = (-jac(0,2)*jac(1,1)+jac(0,1)*jac(1,2))/detjac;//-a02 a11 + a01 a12
1330  jacinv(1,0) = ( jac(1,2)*jac(2,0)-jac(1,0)*jac(2,2))/detjac;//a12 a20 - a10 a22
1331  jacinv(1,1) = (-jac(0,2)*jac(2,0)+jac(0,0)*jac(2,2))/detjac;//-a02 a20 + a00 a22
1332  jacinv(1,2) = ( jac(0,2)*jac(1,0)-jac(0,0)*jac(1,2))/detjac;//a02 a10 - a00 a12
1333  jacinv(2,0) = (-jac(1,1)*jac(2,0)+jac(1,0)*jac(2,1))/detjac;//-a11 a20 + a10 a21
1334  jacinv(2,1) = ( jac(0,1)*jac(2,0)-jac(0,0)*jac(2,1))/detjac;//a01 a20 - a00 a21
1335  jacinv(2,2) = (-jac(0,1)*jac(1,0)+jac(0,0)*jac(1,1))/detjac;//-a01 a10 + a00 a11
1336 
1337 
1338  axes.Zero();
1339  axes(0,0) = 1.0;
1340  axes(1,1) = 1.0;
1341  axes(2,2) = 1.0;
1342 
1343  }
1344  break;
1345 
1346  default:
1347  {
1348  std::cout << " Object with wrong dimensions, unable to compute jacobian matrix. Dimension = " << dim << std::endl;
1349  DebugStop();
1350  }
1351  break;
1352  }
1353 
1354 
1355 }
1356 
1357 void TPZGeoEl::JacobianXYZ(const TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &jac,TPZFMatrix<REAL> &axesXYZ,REAL &detjac,TPZFMatrix<REAL> &jacinv){
1358 
1359  int nrows = gradx.Rows();
1360  int ncols = gradx.Cols();
1361  int dim = ncols;
1362 
1363  switch (dim) {
1364  case 1:
1365  {
1366  jac.Resize(dim,dim);
1367  axesXYZ.Resize(dim,3);
1368  jacinv.Resize(dim,dim);
1369  jac.Zero();
1370 
1371 
1372  if (gradx.GetVal(0,0) != 0.0) { // X oriented
1373 #ifdef PZDEBUG
1374  if (! ( IsZero(gradx.GetVal(1,0)) && IsZero(gradx.GetVal(2,0)) ) ) {
1375  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with x axis. Should be called TPZGeoEl::Jacobian " << std::endl;
1376  DebugStop();
1377  }
1378 #endif
1379  jac(0,0) = gradx.GetVal(0,0);
1380  detjac = gradx.GetVal(0,0);
1381  jacinv(0,0) = 1.0/gradx.GetVal(0,0);
1382 
1383  axesXYZ.Zero();
1384  axesXYZ(0,0) = 1.0;
1385 
1386  if(IsZero(detjac))
1387  {
1388 
1389 #ifdef PZDEBUG
1390  std::stringstream sout;
1391  sout << "Singular Jacobian, 4 determinant of jacobian = " << detjac << std::endl;
1392  LOGPZ_ERROR(logger, sout.str())
1393  DebugStop();
1394 #endif
1395  detjac = ZeroTolerance();
1396  }
1397 
1398  return;
1399  }
1400 
1401 
1402  if (gradx.GetVal(1,0) != 0.0) { // Y oriented
1403 #ifdef PZDEBUG
1404  if (! ( IsZero(gradx.GetVal(0,0)) && IsZero(gradx.GetVal(2,0)) ) ) {
1405  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with y axis. Should be called TPZGeoEl::Jacobian " << std::endl;
1406  DebugStop();
1407  }
1408 #endif
1409  jac(0,0) = gradx.GetVal(1,0);
1410  detjac = gradx.GetVal(1,0);
1411  jacinv(0,0) = 1.0/gradx.GetVal(1,0);
1412 
1413  axesXYZ.Zero();
1414  axesXYZ(0,1) = 1.0;
1415 
1416  if(IsZero(detjac))
1417  {
1418 
1419 #ifdef PZDEBUG
1420  std::stringstream sout;
1421  sout << "Singular Jacobian, 5 determinant of jacobian = " << detjac << std::endl;
1422  LOGPZ_ERROR(logger, sout.str())
1423  DebugStop();
1424 #endif
1425  detjac = ZeroTolerance();
1426  }
1427 
1428  return;
1429  }
1430 
1431  if (gradx.GetVal(2,0) != 0.0) { // Z oriented
1432 #ifdef PZDEBUG
1433  if (! ( IsZero(gradx.GetVal(0,0)) && IsZero(gradx.GetVal(1,0)) ) ) {
1434  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with z axis. Should be called TPZGeoEl::Jacobian " << std::endl;
1435  DebugStop();
1436  }
1437 #endif
1438  jac(0,0) = gradx.GetVal(2,0);
1439  detjac = gradx.GetVal(2,0);
1440  jacinv(0,0) = 1.0/gradx.GetVal(2,0);
1441 
1442  axesXYZ.Zero();
1443  axesXYZ(0,2) = 1.0;
1444 
1445  if(IsZero(detjac))
1446  {
1447 
1448 #ifdef PZDEBUG
1449  std::stringstream sout;
1450  sout << "Singular Jacobian, 6 determinant of jacobian = " << detjac << std::endl;
1451  LOGPZ_ERROR(logger, sout.str())
1452  DebugStop();
1453 #endif
1454  detjac = ZeroTolerance();
1455  }
1456 
1457  return;
1458  }
1459 
1460  }
1461  break;
1462  case 2:
1463  {
1464  jac.Resize(dim,dim);
1465  axesXYZ.Resize(dim,3);
1466  jacinv.Resize(dim,dim);
1467  jac.Zero();
1468 
1469  if (IsZero(gradx.GetVal(2,0)) && IsZero(gradx.GetVal(2,1))) { // XY oriented
1470 #ifdef PZDEBUG
1471  if ( ( IsZero(gradx.GetVal(0,0)) && IsZero(gradx.GetVal(0,1)) ) ) {
1472  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with XY plane. Should be called TPZGeoEl::Jacobian " << std::endl;
1473  DebugStop();
1474  }
1475 #endif
1476 
1477  jac(0,0) = gradx.GetVal(0,0);
1478  jac(1,0) = gradx.GetVal(1,0);
1479  jac(0,1) = gradx.GetVal(0,1);
1480  jac(1,1) = gradx.GetVal(1,1);
1481 
1482  detjac = jac(0,0)*jac(1,1)-jac(1,0)*jac(0,1);
1483  jacinv(0,0) = +jac(1,1)/detjac;
1484  jacinv(1,1) = +jac(0,0)/detjac;
1485  jacinv(0,1) = -jac(0,1)/detjac;
1486  jacinv(1,0) = -jac(1,0)/detjac;
1487 
1488  axesXYZ.Zero();
1489  axesXYZ(0,0) = 1.0;
1490  axesXYZ(1,1) = 1.0;
1491 
1492  if(IsZero(detjac))
1493  {
1494 
1495 #ifdef PZDEBUG
1496  std::stringstream sout;
1497  sout << "Singular Jacobian, 7 determinant of jacobian = " << detjac << std::endl;
1498  LOGPZ_ERROR(logger, sout.str())
1499  DebugStop();
1500 #endif
1501  detjac = ZeroTolerance();
1502  }
1503 
1504  return;
1505  }
1506 
1507  if (IsZero(gradx.GetVal(1,0)) && IsZero(gradx.GetVal(1,1))) { // XZ oriented
1508 #ifdef PZDEBUG
1509  if ( ( IsZero(gradx.GetVal(0,0)) && IsZero(gradx.GetVal(0,1)) ) ) {
1510  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with XZ plane. Should be called TPZGeoEl::Jacobian " << std::endl;
1511  DebugStop();
1512  }
1513 #endif
1514 
1515  jac(0,0) = gradx.GetVal(0,0);
1516  jac(1,0) = gradx.GetVal(2,0);
1517  jac(0,1) = gradx.GetVal(0,1);
1518  jac(1,1) = gradx.GetVal(2,1);
1519 
1520  detjac = jac(0,0)*jac(1,1)-jac(1,0)*jac(0,1);
1521  jacinv(0,0) = +jac(1,1)/detjac;
1522  jacinv(1,1) = +jac(0,0)/detjac;
1523  jacinv(0,1) = -jac(0,1)/detjac;
1524  jacinv(1,0) = -jac(1,0)/detjac;
1525 
1526  axesXYZ.Zero();
1527  axesXYZ(0,0) = 1.0;
1528  axesXYZ(1,2) = 1.0;
1529 
1530  if(IsZero(detjac))
1531  {
1532 
1533 #ifdef PZDEBUG
1534  std::stringstream sout;
1535  sout << "Singular Jacobian, 8 determinant of jacobian = " << detjac << std::endl;
1536  LOGPZ_ERROR(logger, sout.str())
1537  DebugStop();
1538 #endif
1539  detjac = ZeroTolerance();
1540  }
1541 
1542  return;
1543  }
1544 
1545  if (IsZero(gradx.GetVal(0,0)) && IsZero(gradx.GetVal(0,1))) { // YZ oriented
1546 #ifdef PZDEBUG
1547  if ( ( IsZero(gradx.GetVal(1,0)) && IsZero(gradx.GetVal(1,1)) ) ) {
1548  std::cout << "TPZGeoEl::JacobianXYZ -> Geometric mesh is not oriented with YZ plane. Should be called TPZGeoEl::Jacobian " << std::endl;
1549  DebugStop();
1550  }
1551 #endif
1552 
1553  jac(0,0) = gradx.GetVal(1,0);
1554  jac(1,0) = gradx.GetVal(2,0);
1555  jac(0,1) = gradx.GetVal(1,1);
1556  jac(1,1) = gradx.GetVal(2,1);
1557 
1558  detjac = jac(0,0)*jac(1,1)-jac(1,0)*jac(0,1);
1559  jacinv(0,0) = +jac(1,1)/detjac;
1560  jacinv(1,1) = +jac(0,0)/detjac;
1561  jacinv(0,1) = -jac(0,1)/detjac;
1562  jacinv(1,0) = -jac(1,0)/detjac;
1563 
1564  axesXYZ.Zero();
1565  axesXYZ(0,1) = 1.0;
1566  axesXYZ(1,2) = 1.0;
1567 
1568  if(IsZero(detjac))
1569  {
1570 
1571 #ifdef PZDEBUG
1572  std::stringstream sout;
1573  sout << "Singular Jacobian, 9 determinant of jacobian = " << detjac << std::endl;
1574  LOGPZ_ERROR(logger, sout.str())
1575  DebugStop();
1576 #endif
1577  detjac = ZeroTolerance();
1578  }
1579 
1580  return;
1581  }
1582 
1583 
1584  }
1585  break;
1586  case 3:
1587  {
1588  jac.Resize(dim,dim);
1589  axesXYZ.Resize(dim,3);
1590  jacinv.Resize(dim,dim);
1591  jac.Zero();
1592 
1593  for (int i = 0; i < nrows; i++) {
1594  jac(i,0) = gradx.GetVal(i,0);
1595  jac(i,1) = gradx.GetVal(i,1);
1596  jac(i,2) = gradx.GetVal(i,2);
1597  }
1598 
1599  detjac = -jac(0,2)*jac(1,1)*jac(2,0);//-a02 a11 a20
1600  detjac += jac(0,1)*jac(1,2)*jac(2,0);//+ a01 a12 a20
1601  detjac += jac(0,2)*jac(1,0)*jac(2,1);//+ a02 a10 a21
1602  detjac -= jac(0,0)*jac(1,2)*jac(2,1);//- a00 a12 a21
1603  detjac -= jac(0,1)*jac(1,0)*jac(2,2);//- a01 a10 a22
1604  detjac += jac(0,0)*jac(1,1)*jac(2,2);//+ a00 a11 a22
1605 
1606  jacinv(0,0) = (-jac(1,2)*jac(2,1)+jac(1,1)*jac(2,2))/detjac;//-a12 a21 + a11 a22
1607  jacinv(0,1) = ( jac(0,2)*jac(2,1)-jac(0,1)*jac(2,2))/detjac;//a02 a21 - a01 a22
1608  jacinv(0,2) = (-jac(0,2)*jac(1,1)+jac(0,1)*jac(1,2))/detjac;//-a02 a11 + a01 a12
1609  jacinv(1,0) = ( jac(1,2)*jac(2,0)-jac(1,0)*jac(2,2))/detjac;//a12 a20 - a10 a22
1610  jacinv(1,1) = (-jac(0,2)*jac(2,0)+jac(0,0)*jac(2,2))/detjac;//-a02 a20 + a00 a22
1611  jacinv(1,2) = ( jac(0,2)*jac(1,0)-jac(0,0)*jac(1,2))/detjac;//a02 a10 - a00 a12
1612  jacinv(2,0) = (-jac(1,1)*jac(2,0)+jac(1,0)*jac(2,1))/detjac;//-a11 a20 + a10 a21
1613  jacinv(2,1) = ( jac(0,1)*jac(2,0)-jac(0,0)*jac(2,1))/detjac;//a01 a20 - a00 a21
1614  jacinv(2,2) = (-jac(0,1)*jac(1,0)+jac(0,0)*jac(1,1))/detjac;//-a01 a10 + a00 a11
1615 
1616  if(IsZero(detjac))
1617  {
1618 
1619 #ifdef PZDEBUG
1620  std::stringstream sout;
1621  sout << "Singular Jacobian, 10 determinant of jacobian = " << detjac << std::endl;
1622  LOGPZ_ERROR(logger, sout.str())
1623  DebugStop();
1624 #endif
1625  detjac = ZeroTolerance();
1626  }
1627 
1628  axesXYZ.Zero();
1629  axesXYZ(0,0) = 1.0;
1630  axesXYZ(1,1) = 1.0;
1631  axesXYZ(2,2) = 1.0;
1632 
1633  }
1634  break;
1635 
1636  default:
1637  {
1638  std::cout << " Object with wrong dimensions, unable to compute jacobian matrix. Dimension = " << dim << std::endl;
1639  DebugStop();
1640  }
1641  break;
1642  }
1643 
1644 }
1645 
1647 void TPZGeoEl::SetRefPattern(TPZAutoPointer<TPZRefPattern> ){
1648  PZError << "TPZGeoEl::SetRefPattern ERROR : Should not be called in TPZGeoEl" << endl;
1649  DebugStop();
1650 }
1651 
1652 void TPZGeoEl::Read(TPZStream &buf, void *context) { //ok
1653  fMesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::GetInstance(&buf));
1654  buf.Read(&fId,1);
1655  buf.Read(&fMatId,1);
1656  fReference = dynamic_cast<TPZCompEl *>(TPZPersistenceManager::GetInstance(&buf));
1657  buf.Read(&fFatherIndex,1);
1658  buf.Read(&fIndex,1);
1659  gGlobalAxes.Read(buf, 0);
1660  buf.Read(&fNumInterfaces,1);
1661 }
1662 
1663 void TPZGeoEl::Write(TPZStream &buf, int withclassid) const { //ok
1664  TPZPersistenceManager::WritePointer(fMesh, &buf);
1665  buf.Write(&fId, 1);
1666  buf.Write(&fMatId, 1);
1667  TPZPersistenceManager::WritePointer(fReference, &buf);
1668  buf.Write(&fFatherIndex, 1);
1669  buf.Write(&fIndex, 1);
1670  gGlobalAxes.Write(buf, 0);
1671  buf.Write(&fNumInterfaces, 1);
1672 }
1673 
1674 TPZGeoEl::TPZGeoEl(TPZGeoMesh & DestMesh, const TPZGeoEl &cp):TPZSavable(cp){
1675  this->fMesh = &DestMesh;
1676  this->fId = cp.fId;
1677  this->fMatId = cp.fMatId;
1678  this->fReference = cp.fReference;
1679  this->fFatherIndex = cp.fFatherIndex;
1680  this->fIndex = cp.fIndex;
1681  this->fMesh->ElementVec()[this->fIndex] = this;
1682  this->fNumInterfaces = 0;
1683 }
1684 
1685 TPZGeoEl::TPZGeoEl(TPZGeoMesh & DestMesh, const TPZGeoEl &cp, std::map<int64_t,int64_t> &org2clnMap):TPZSavable(cp){
1686  this->fMesh = &DestMesh;
1687  this->fId = cp.fId;
1688  this->fMatId = cp.fMatId;
1689  this->fReference = 0;
1690  if ( cp.fFatherIndex == -1) this->fFatherIndex = -1;
1691  else if (org2clnMap.find(cp.fFatherIndex) == org2clnMap.end())
1692  {
1693  std::stringstream sout;
1694  sout << "ERROR in - " << __PRETTY_FUNCTION__
1695  << " original father element index: " << cp.fIndex << " is not mapped!";
1696  LOGPZ_ERROR (logger,sout.str().c_str());
1697  DebugStop();
1698  }
1699  else this->fFatherIndex = org2clnMap[cp.fFatherIndex];
1700 
1701  if (org2clnMap.find(cp.fIndex) == org2clnMap.end())
1702  {
1703  std::stringstream sout;
1704  sout << "ERROR in - " << __PRETTY_FUNCTION__
1705  << " original element index: " << cp.fIndex << " is not mapped!";
1706  LOGPZ_ERROR (logger,sout.str().c_str());
1707  DebugStop();
1708  }
1709 
1710  this->fIndex = org2clnMap[cp.fIndex];
1711  this->fMesh->ElementVec()[this->fIndex] = this;
1712  this->fNumInterfaces = 0;
1713 }
1714 
1715 // return the refinement pattern associated with the element
1716 TPZAutoPointer<TPZRefPattern> TPZGeoEl::GetRefPattern() const
1717 {
1718  TPZAutoPointer<TPZRefPattern> result;
1719  return result;
1720 }
1721 
1722 bool TPZGeoEl::VerifyNodeCoordinates(REAL tol){
1723  const int nnodes = this->NCornerNodes();
1724  TPZManVector<REAL,3> qsi(this->Dimension());
1725  TPZManVector<REAL,3> MappedX(3), NodeX(3);
1726  for(int inode = 0; inode < nnodes; inode++){
1727  for(int dim = 0; dim < 3; dim++){
1728  NodeX[dim] = this->NodePtr(inode)->Coord(dim);
1729  }//dim
1730  this->CenterPoint(inode,qsi);
1731  this->X(qsi,MappedX);
1732  REAL error = 0.;
1733  for(int dim = 0; dim < 3; dim++){
1734  error += (NodeX[dim]-MappedX[dim])*(NodeX[dim]-MappedX[dim]);
1735  }//dim
1736  error = sqrt(error);
1737  if(error > tol || !(error==error)){
1738  std::stringstream mess;
1739  mess << "FATAL ERROR AT " << __PRETTY_FUNCTION__ << " - Node coordinate differs from mapped node.\n";
1740  this->Print(mess);
1741  PZError << mess.str() << "\n";
1742 #ifdef LOG4CXX
1743  LOGPZ_ERROR(logger,mess.str().c_str());
1744 #endif
1745  DebugStop();
1746  return false;
1747  }
1748  }//for i
1749  return true;
1750 }//method
1751 
1752 
1753 #include "tpzgeoelrefpattern.h"
1754 #include "TPZGeoCube.h"
1755 #include "pzshapecube.h"
1756 #include "TPZRefCube.h"
1757 #include "pzshapelinear.h"
1758 #include "TPZGeoLinear.h"
1759 #include "TPZRefLinear.h"
1760 #include "pzrefquad.h"
1761 #include "pzshapequad.h"
1762 #include "pzgeoquad.h"
1763 #include "pzshapetriang.h"
1764 #include "pzreftriangle.h"
1765 #include "pzgeotriangle.h"
1766 #include "pzshapeprism.h"
1767 #include "pzrefprism.h"
1768 #include "pzgeoprism.h"
1769 #include "pzshapetetra.h"
1770 #include "pzreftetrahedra.h"
1771 #include "pzgeotetrahedra.h"
1772 #include "pzshapepiram.h"
1773 #include "pzrefpyram.h"
1774 #include "pzgeopyramid.h"
1775 #include "pzgeopoint.h"
1776 #include "pzrefpoint.h"
1777 #include "pzshapepoint.h"
1778 
1779 using namespace pzgeom;
1780 using namespace pzrefine;
1781 using namespace pzshape;
1782 
1783 
1785 int ConjugateSide(TPZGeoEl *gel, int side, TPZStack<int> &allsides);
1787 void NormalVector(TPZGeoElSide &LC, TPZGeoElSide &LS, TPZVec<REAL> &normal);
1789 void Normalize(TPZVec<REAL> &normlow, TPZVec<REAL> &normal);
1790 
1791 //void TPZGeoEl::ComputeNormals(TPZMatrix<REAL> &normals)
1792 //{
1793 //
1794 //#ifdef LOG4CXX
1795 // {
1796 // std::stringstream sout;
1797 // sout<< "Metodo Compute normal \n";
1798 //
1799 // LOGPZ_DEBUG(logger,sout.str())
1800 // }
1801 //#endif
1802 //
1803 //
1804 // int numbernormals = 0;
1805 // int dimension = Dimension();
1806 // int is;
1807 // int nsides = NSides();
1808 // for(is=0; is<nsides; is++)
1809 // {
1810 // TPZStack<int> lowdim;
1811 // if(SideDimension(is) == dimension-1)
1812 // {
1813 // //TPZStack<int> lowdim;
1814 // LowerDimensionSides(is,lowdim);
1815 // numbernormals += lowdim.NElements()+1;
1816 // }
1817 //
1818 // }
1819 // normals.Redim(3, numbernormals);
1820 // int counter = 0;
1821 // for(is=0; is<nsides; is++)
1822 // {
1823 // if(SideDimension(is) == dimension-1)
1824 //
1825 // {
1826 //
1827 //#ifdef LOG4CXX
1828 // {
1829 // std::stringstream sout;
1830 // sout<< "Side "<<is<<std::endl;
1831 //
1832 // LOGPZ_DEBUG(logger,sout.str())
1833 // }
1834 //#endif
1835 //
1836 //
1837 // TPZStack<int> lowdim;
1838 // LowerDimensionSides(is,lowdim);
1839 // lowdim.Push(is);
1840 //#ifdef LOG4CXX
1841 // {
1842 // std::stringstream sout;
1843 // sout<< "LowerDimensionSides "<<lowdim<<std::endl;
1844 //
1845 // LOGPZ_DEBUG(logger,sout.str())
1846 // }
1847 //#endif
1848 //
1849 // int nlowdim = lowdim.NElements();
1850 // int lowis;
1851 // for(lowis=0; lowis < nlowdim; lowis++)
1852 // {
1853 // int conj_side = ConjugateSide(this,lowdim[lowis],lowdim);
1854 // TPZGeoElSide LC(this,conj_side);
1855 // TPZGeoElSide LS(this,lowdim[lowis]);
1856 //#ifdef LOG4CXX
1857 // {
1858 // std::stringstream sout;
1859 // sout<< "Side "<<is << " Conjugate Side "<< conj_side<< " Ls side "<< lowdim[lowis]<<std::endl;
1860 //
1861 // LOGPZ_DEBUG(logger,sout.str())
1862 // }
1863 //#endif
1864 // TPZManVector<REAL> normal(3,0.);
1865 // NormalVector(LC,LS,normal);
1866 //
1867 //#ifdef LOG4CXX
1868 // {
1869 // std::stringstream sout;
1870 // sout<< "Vetores do NormalVector "<<normal<<std::endl;
1871 //
1872 // LOGPZ_DEBUG(logger,sout.str())
1873 // }
1874 //#endif
1875 //
1876 // int d;
1877 // for(d=0; d<3; d++) normals(d,counter) = normal[d];
1878 // counter++;
1879 // }
1880 //
1881 //#ifdef LOG4CXX
1882 // {
1883 // std::stringstream sout;
1884 // sout<< "A partir daqui sera um processo de normalizacao dos vetores normals"<<std::endl;
1885 // normals.Print("vetor normals",sout);
1886 //
1887 // LOGPZ_DEBUG(logger,sout.str())
1888 // }
1889 //#endif
1890 // TPZManVector<REAL,3> normal(3,0.);
1891 // int d;
1892 // for(d=0; d<3; d++) normal[d] = normals(d,counter-1);
1893 // for(lowis = counter-nlowdim; lowis < counter; lowis++)
1894 // {
1895 // TPZManVector<REAL,3> normlow(3,0.);
1896 // for(d=0; d<3; d++) normlow[d] = normals(d,lowis);
1897 // Normalize(normlow,normal);
1898 // for(d=0; d<3; d++) normals(d,lowis) = normlow[d];
1899 // }
1900 // }
1901 // }
1902 //
1903 //#ifdef LOG4CXX
1904 // {
1905 // std::stringstream sout;
1906 // sout<< "Vetores Normais normalizados "<<normals<<std::endl;
1907 //
1908 // LOGPZ_DEBUG(logger,sout.str())
1909 // }
1910 //#endif
1911 //}
1912 
1913 void TPZGeoEl::SetNeighbourForBlending(int side){
1914  if( !this->IsGeoBlendEl() ) return;
1915 
1916  TPZGeoElSide ElemSide(this,side);
1917  TPZGeoElSide NextSide(this,side);
1918 
1919  //trying to get a self non-linear neighbour
1920  TPZGeoElSide NextSideNeighbour = NextSide.Neighbour();
1921  while(NextSideNeighbour.Element() != ElemSide.Element())
1922  {
1923  if(NextSideNeighbour.Exists() && !NextSideNeighbour.Element()->IsLinearMapping() && !NextSideNeighbour.Element()->IsGeoBlendEl())
1924  {
1925  if(NextSideNeighbour.IsRelative(ElemSide) == false){
1926  if(NextSideNeighbour.Element()->IsGeoElMapped() == false){
1927  TPZGeoElSide NeighSide = NextSideNeighbour;
1928  TPZTransform<> NeighTransf(NeighSide.Dimension(),NeighSide.Dimension());
1929  ElemSide.SideTransform3(NeighSide,NeighTransf);
1930  this->SetNeighbourInfo(side,NeighSide,NeighTransf);
1931  return;
1932  }// !TPZGeoElMapped
1933  }//if IsRelative == false
1934  }
1935  NextSide = NextSide.Neighbour();
1936  NextSideNeighbour = NextSide.Neighbour();
1937  }
1938  NextSide = TPZGeoElSide(this,side);
1939  NextSideNeighbour = NextSide.Neighbour();
1940  //now TPZGeoElMapped are accepted
1941  while(NextSideNeighbour.Element() != ElemSide.Element())
1942  {
1943  if(NextSideNeighbour.Exists() && !NextSideNeighbour.Element()->IsLinearMapping() && !NextSideNeighbour.Element()->IsGeoBlendEl())
1944  {
1945  if(NextSideNeighbour.IsRelative(ElemSide) == false){
1946  TPZGeoElSide NeighSide = NextSideNeighbour;
1947  TPZTransform<> NeighTransf(NeighSide.Dimension(),NeighSide.Dimension());
1948  ElemSide.SideTransform3(NeighSide,NeighTransf);
1949  this->SetNeighbourInfo(side,NeighSide,NeighTransf);
1950  return;
1951  }//if IsRelative == false
1952  }
1953  NextSide = NextSide.Neighbour();
1954  NextSideNeighbour = NextSide.Neighbour();
1955  }
1956 
1957 }//method
1958 
1959 void TPZGeoEl::BuildBlendConnectivity(){
1960  if( !this->IsGeoBlendEl() ) return;
1961  const int nsides = this->NSides();
1962  for(int byside = this->NNodes(); byside < nsides; byside++)
1963  {
1964  this->SetNeighbourForBlending(byside);
1965  }//for byside
1966 }//method
1967 
1968 // Projection of the point to the side
1969 // Compute the projection of the point within the interior of the element to the side of the element
1970 TPZTransform<> TPZGeoEl::Projection(int side)
1971 {
1972  TPZTransform<> tr = SideToSideTransform(NSides()-1,side);
1973  TPZTransform<> tr2 = SideToSideTransform(side,NSides()-1);
1974  TPZTransform<> tr3 = tr2.Multiply(tr);
1975  // std::cout << "interior to side " << side << " trans " << tr << std::endl;
1976  // std::cout << "side to interior " << side << " trans " << tr2 << std::endl;
1977  // std::cout << "side to side " << side << " trans " << tr3 << std::endl;
1978  return tr3;
1979 }
1980 
1981 TPZGeoEl * TPZGeoEl::EldestAncestor() const
1982 {
1983  TPZGeoEl *father = this->Father();
1984  if(!father)
1985  {
1986  return NULL;
1987  }
1988  TPZGeoEl *nextfather = NULL;
1989  if(father) nextfather = father->Father();
1990  while(nextfather)
1991  {
1992  father = nextfather;
1993  nextfather = father->Father();
1994  }
1995 
1996  return father;
1997 }
1998 
1999 void NormalVector(TPZGeoElSide &LC, TPZGeoElSide &LS, TPZVec<REAL> &normal)
2000 {
2001  //TPZGeoEl *gel = LC.Element();
2002  // take the centerpoint of LC and the centerpoint of LS
2003  TPZManVector<REAL,3> LCCenter(3,0.), LSCenter(3,0.);
2004  TPZManVector<REAL,3> XLC(3,0.),XLS(3,0.);
2005  LC.CenterPoint(LCCenter);
2006  LC.X(LCCenter,XLC);
2007 
2008  LS.CenterPoint(LSCenter);
2009  LS.X(LSCenter,XLS);
2010  TPZManVector<REAL,3> dir(3,0.);
2011 
2012 #ifdef LOG4CXX
2013  if(logger->isDebugEnabled())
2014  {
2015  std::stringstream sout;
2016  sout<< "Centro de LC "<<XLC<< " Centro de LS "<<XLS<<std::endl;
2017  LOGPZ_DEBUG(logger,sout.str())
2018  }
2019 #endif
2020 
2021 
2022  // The normal vector needs to be in the plane of LC and perpendicular to LS
2023  // A starting vector is the direction of one center to the next
2024  int i;
2025  for(i=0; i<3; i++) dir[i] = XLS[i]-XLC[i];
2026 
2027 
2028 
2029  TPZFNMatrix<10> jacobian,axes,jacinv;
2030  REAL detjac;
2031  LS.Jacobian(LSCenter,jacobian,axes,detjac,jacinv);
2032 
2033 #ifdef LOG4CXX
2034  if(logger->isDebugEnabled())
2035  {
2036  std::stringstream sout;
2037  sout<< "vetor axes "<<axes<<std::endl;
2038  sout<< "vetor dir "<<dir<<std::endl;
2039  LOGPZ_DEBUG(logger,sout.str())
2040  }
2041 #endif
2042 
2043  TPZFNMatrix<20> axtrans(axes.Cols(),axes.Rows()+1,0.);
2044  int j;
2045  for(i=0; i<3; i++) for(j=0; j<axes.Rows(); j++)
2046  {
2047  axtrans(i,j) = axes(j,i);
2048  }
2049  int lastcol = j;
2050  for(i=0; i<3; i++)
2051  {
2052  axtrans(i,lastcol) = dir[i];
2053  }
2054 
2055 
2056 #ifdef LOG4CXX
2057  if(logger->isDebugEnabled())
2058  {
2059  std::stringstream sout;
2060  sout << "axtrans = ";
2061  axtrans.Print("axtrans matrix",sout);
2062  LOGPZ_DEBUG(logger,sout.str())
2063  }
2064 #endif
2065 
2066 
2067  // Then orthogonalize the vectors of the LS side with this vector
2068  TPZFNMatrix<20> transf, ortho;
2069  axtrans.GramSchmidt(ortho,transf);
2070 
2071 
2072 #ifdef LOG4CXX
2073  if(logger->isDebugEnabled())
2074  {
2075  std::stringstream sout;
2076  sout << "Vetr Apos GramSchmidt = ";
2077  ortho.Print("Orthog matrix",sout);
2078  transf.Print("TransfToOrthog matrix",sout);
2079  LOGPZ_DEBUG(logger,sout.str())
2080  }
2081 #endif
2082 
2083 
2084  normal.Resize(3);
2085  for(i=0; i<3; i++)
2086  {
2087  normal[i] = ortho(i,lastcol);
2088  }
2089 
2090 
2091 }
2092 
2093 void Normalize(TPZVec<REAL> &normlow, TPZVec<REAL> &normal)
2094 {
2095 #ifdef LOG4CXX
2096  if(logger->isDebugEnabled())
2097  {
2098  std::stringstream sout;
2099  sout<< "Normalize..u.n=1 "<<std::endl;
2100  sout<< "normlow "<<normlow<< " normal "<<normal<<std::endl;
2101  LOGPZ_DEBUG(logger,sout.str())
2102  }
2103 #endif
2104 
2105 
2106  REAL inner = 0.;
2107  int i;
2108  for(i=0; i<normlow.NElements(); i++)
2109  {
2110  inner += normlow[i]*normal[i];
2111  }
2112 
2113  if(inner==0)
2114  {
2115  LOGPZ_ERROR(logger,"Inner product zero to normalize vector")
2116  DebugStop();
2117  }
2118 
2119 
2120 
2121  for(i=0; i<normlow.NElements(); i++)
2122  {
2123  normlow[i] /= inner;
2124  }
2125 }
2126 
2127 //void TPZGeoEl::ComputeNormals(int side, TPZFMatrix<REAL> &normals, TPZVec<int> &vectorsides)
2128 //{
2129 // int numbernormals = 0;
2130 // int dimension = Dimension();
2131 // int sidedimension = SideDimension(side);
2132 // TPZStack<int> lowdim;
2133 // LowerDimensionSides(side,lowdim);
2134 // lowdim.Push(side);
2135 // // the normals corresponding to the internal shape functions
2136 // // Compute the number of normals we need to compute
2137 // int nsides = NSides();
2138 // if(sidedimension == dimension-1)
2139 // {
2140 // numbernormals = lowdim.NElements();
2141 // }
2142 // else if(sidedimension == dimension)
2143 // {
2144 // numbernormals = nsides*2;
2145 // }
2146 // else
2147 // {
2148 // numbernormals = 0;
2149 // }
2150 // normals.Redim(3, numbernormals);
2151 // vectorsides.Resize(numbernormals);
2152 // vectorsides.Fill(0);
2153 // if(!numbernormals)
2154 // {
2155 // return;
2156 // }
2157 // int is = side;
2158 // if(sidedimension == dimension-1)
2159 // {
2160 // int nlowdim = lowdim.NElements();
2161 // int lowis;
2162 // // work from lowest dimension sides upward
2163 // for(lowis=0; lowis < nlowdim; lowis++)
2164 // {
2165 // // find a side which is not contained in the currently analysed side
2166 // // whose dimension is one higher than the dimension of lowdim[lowis]
2167 // int conj_side = ConjugateSide(this,lowdim[lowis],lowdim);
2168 // // the normal vector is in the alignment of the conjugate side
2169 // TPZGeoElSide LC(this,conj_side);
2170 // TPZGeoElSide LS(this,lowdim[lowis]);
2171 // TPZManVector<REAL> normal(3,0.);
2172 // // the normal vector goes from the center of the conjugate side to
2173 // // the center of the LS side
2174 // NormalVector(LC,LS,normal);
2175 // int d;
2176 // for(d=0; d<3; d++) normals(d,lowis) = normal[d];
2177 // vectorsides[lowis] = lowdim[lowis];
2178 // }
2179 // // Why mormalize a vector which has just been computed?
2180 // TPZManVector<REAL> normal(3,0.);
2181 // int d;
2182 // for(d=0; d<3; d++) normal[d] = normals(d,nlowdim-1);
2183 // for(lowis = 0; lowis < nlowdim; lowis++)
2184 // {
2185 // TPZManVector<REAL,3> normlow(3,0.);
2186 // for(d=0; d<3; d++) normlow[d] = normals(d,lowis);
2187 // Normalize(normlow,normal);
2188 // for(d=0; d<3; d++) normals(d,lowis) = normlow[d];
2189 // }
2190 // TPZManVector<int,9> sidepermutationgather(nlowdim);
2192 // for(int i=0; i< nlowdim; i++) sidepermutationgather[i] = i;
2193 // TPZFNMatrix<12> sidenormals(3,nlowdim);
2194 // TPZManVector<int> localvecsides(nlowdim);
2195 // // compute whether the side is from this element to the next or contrary
2196 // int sideorient = 1;//NormalOrientation(side);
2197 // int i;
2198 // for(i=0; i<nlowdim; i++)
2199 // {
2200 // for(d=0; d<3; d++)
2201 // {
2202 // sidenormals(d,i) = normals(d,sidepermutationgather[i]);
2203 // }
2204 // localvecsides[i] = vectorsides[sidepermutationgather[i]];
2205 // }
2206 // for(i=0; i<nlowdim; i++)
2207 // {
2208 // for(d=0; d<3; d++)
2209 // {
2210 // normals(d,i) = sidenormals(d,i)*sideorient;
2211 // }
2212 // vectorsides[i] = localvecsides[i];
2213 // }
2214 // }
2215 // else if(sidedimension == dimension)
2216 // {
2217 // int counter = 0;
2218 // for(is=0; is<nsides; is++)
2219 // {
2220 // if(SideDimension(is) > 0)
2221 // {
2222 // TPZManVector<REAL,3> Center(3,0.);
2223 // TPZManVector<REAL,3> X;
2224 // TPZFNMatrix<10> jacobian, jacinv, axes;
2225 // REAL detjac;
2226 // TPZGeoElSide gelside(this,is);
2227 // gelside.CenterPoint(Center);
2228 // gelside.Jacobian(Center,jacobian,axes,detjac,jacinv);
2229 // int d,s;
2230 // dimension = SideDimension(is);
2231 // for(d=0; d<dimension; d++)
2232 // {
2233 // for(s=0; s<3; s++)
2234 // {
2235 // normals(s,counter) = axes(d,s);
2236 // }
2237 // vectorsides[counter] = is;
2238 // counter++;
2239 // }
2240 // }
2241 // }
2242 // vectorsides.Resize(counter);
2243 // normals.Resize(3,counter);
2244 // }
2245 //}
2246 
2247 
2248 //void TPZGeoEl::ComputeNormalsDG(int side, TPZVec<REAL> &pt, TPZFMatrix<REAL> &normals, TPZVec<int> &vectorsides)
2249 //{
2250 // if (SideDimension(side) >= Dimension()-1) {
2251 // Directions(side, pt, normals, vectorsides);
2252 // }
2253 // if (SideDimension(side) == Dimension()-1)
2254 // {
2255 // // we need to permute the normals and associated sides
2256 // TPZGeoElSide thisside(this,side);
2257 // int nlowdim = thisside.NSides();
2258 //
2259 // TPZManVector<int,9> sidepermutationgather(nlowdim);
2260 // HDivPermutation(side,sidepermutationgather);
2261 //#ifdef LOG4CXX
2262 // if(logger->isDebugEnabled()){
2263 // std::stringstream sout;
2264 // sout << "Permutation for side " << side << " is " << sidepermutationgather;
2265 // LOGPZ_DEBUG(logger, sout.str())
2266 // }
2267 //#endif
2268 // TPZFNMatrix<12> sidenormals(3,nlowdim);
2269 // TPZManVector<int> localvecsides(nlowdim);
2270 // // compute whether the side is from this element to the next or contrary
2271 // int sideorient = NormalOrientation(side);
2272 // int i;
2273 // for(i=0; i<nlowdim; i++)
2274 // {
2275 // for(int d=0; d<3; d++)
2276 // {
2277 // sidenormals(d,i) = normals(d,sidepermutationgather[i]);
2278 // }
2279 // localvecsides[i] = vectorsides[sidepermutationgather[i]];
2280 // }
2281 // for(i=0; i<nlowdim; i++)
2282 // {
2283 // for(int d=0; d<3; d++)
2284 // {
2285 // normals(d,i) = sidenormals(d,i)*sideorient;
2286 // }
2287 // vectorsides[i] = localvecsides[i];
2288 // }
2289 //
2290 // }
2291 //
2292 //}
2293 
2294 //void TPZGeoEl::ComputeNormals(TPZFMatrix<REAL> &normals, TPZVec<int> &vectorsides)
2295 //{
2296 // int numbernormals = 0;
2297 // int dimension = Dimension();
2298 // // the normals corresponding to the internal shape functions
2299 // int is;
2300 // // Compute the number of normals we need to compute
2301 // int nsides = NSides();
2302 // numbernormals = nsides*dimension; // @omar:: why two???
2303 // normals.Redim(3, numbernormals);
2304 // vectorsides.Resize(numbernormals);
2305 // vectorsides.Fill(0);
2306 // int counter = 0;
2307 // // effectively compute the normals
2308 // for(is=0; is<nsides; is++)
2309 // {
2310 // TPZFNMatrix<100> sidenormals;
2311 // TPZManVector<int> sidevectors;
2312 // ComputeNormals(is,sidenormals,sidevectors);
2313 // int numnormals = sidevectors.NElements();
2314 // int in;
2315 // for(in=0; in<numnormals; in++)
2316 // {
2317 // int d;
2318 // for(d=0; d<3; d++)
2319 // {
2320 // normals(d,counter) = sidenormals(d,in);
2321 // }
2322 // vectorsides[counter] = sidevectors[in];
2323 // counter++;
2324 // }
2325 // }
2326 //#ifdef PZDEBUG
2327 // if(counter != numbernormals)
2328 // {
2329 // DebugStop();
2330 // }
2331 //#endif
2332 //}
2333 
2334 //void TPZGeoEl::ComputeNormalsDG(TPZVec<REAL> &pt, TPZFMatrix<REAL> &normals, TPZVec<int> &vectorsides)
2335 //{
2336 // int numbernormals = 0;
2337 // // int dimension = Dimension();
2338 // // the normals corresponding to the internal shape functions
2339 // int is;
2340 // // Compute the number of normals we need to compute
2341 // int nsides = NSides();
2342 // numbernormals = nsides*Dimension();
2343 // normals.Redim(3, numbernormals);
2344 // vectorsides.Resize(numbernormals);
2345 // vectorsides.Fill(0);
2346 // int counter = 0;
2347 // // effectively compute the normals
2348 // for(is=0; is<nsides; is++)
2349 // {
2350 // TPZFNMatrix<100> sidenormals;
2351 // TPZManVector<int> sidevectors;
2352 // ComputeNormalsDG(is,pt, sidenormals,sidevectors);
2353 // int numnormals = sidevectors.NElements();
2354 // int in;
2355 // for(in=0; in<numnormals; in++)
2356 // {
2357 // int d;
2358 // for(d=0; d<3; d++)
2359 // {
2360 // normals(d,counter) = sidenormals(d,in);
2361 // }
2362 // vectorsides[counter] = sidevectors[in];
2363 // counter++;
2364 // }
2365 // }
2366 //}
2367 
2368 
2369 // Determine the orientation of the normal vector comparing the ids of the neighbouring elements
2370 int TPZGeoEl::NormalOrientation(int side)
2371 {
2372  int dimel = Dimension();
2373  int dimside = SideDimension(side);
2374  if(dimside != dimel-1)
2375  {
2376  LOGPZ_ERROR(logger,"NormalOrientation called with wrong side")
2377  return 0;
2378  }
2379 #ifdef LOG4CXX
2380  if (loggerorient->isDebugEnabled())
2381  {
2382  std::stringstream sout;
2383  TPZCompEl *cel = Reference();
2384  if (cel) {
2385  sout << "Comp Element index " << cel->Index() << std::endl;
2386  }
2387  sout << "Element index " << this->Index() << " Id = " << this->Id() << " side = " << side;
2388  LOGPZ_DEBUG(loggerorient, sout.str())
2389  }
2390 #endif
2391  TPZGeoElSide thisside(this,side);
2392 
2393  TPZGeoElSide fatherside = thisside.Father2();
2394  while(fatherside.Exists() && fatherside.Dimension() == dimside)//a segunda condicional eu inclui agora
2395  {
2396  thisside = fatherside;
2397  fatherside = fatherside.Father2();
2398  }
2399 
2400 #ifdef LOG4CXX
2401  if (loggerorient->isDebugEnabled())
2402  {
2403  std::stringstream sout;
2404  sout << "largest father index " << thisside.Element()->Index() << " id = " << thisside.Element()->Id() << std::endl;
2405  LOGPZ_DEBUG(loggerorient, sout.str())
2406  }
2407 #endif
2408 
2409 
2410  TPZGeoElSide neighbour(thisside.Neighbour());
2411  // case there is no neighbour (should put a debugstop here)
2412  if(!neighbour.Exists() )
2413  {
2414  return 1;
2415  }
2416  fatherside = neighbour.Father2();
2417  // look for a neighbour of equal dimension
2418  while ((neighbour.Element()->Dimension() != Dimension() && neighbour != thisside) || (fatherside && fatherside.Dimension() == dimside)) {
2419  neighbour = neighbour.Neighbour();
2420  fatherside = neighbour.Father2();
2421  }
2422  if (neighbour == thisside) {
2423  return 1;
2424  }
2425 
2426 #ifdef PZDEBUG
2427  if(!thisside.NeighbourExists(neighbour))//inclui agora esta verificacao
2428  {
2429  std::stringstream sout;
2430 #ifdef LOG4CXX
2431  LOGPZ_ERROR(logger,sout.str().c_str());
2432 #endif
2433  DebugStop();
2434 
2435  }
2436 #endif
2437 #ifdef LOG4CXX
2438  if (loggerorient->isDebugEnabled())
2439  {
2440  std::stringstream sout;
2441  sout << "Element index " << Index() << std::endl;
2442  sout << "neighbour index " << neighbour.Element()->Index() << " id = " << neighbour.Element()->Id() << " side " << neighbour.Side() << std::endl;
2443  sout << "thisside index " << thisside.Element()->Index() << " id = " << thisside.Element()->Id() << " side " << thisside.Side() << std::endl;
2444  if(thisside.Element()->Id() < neighbour.Element()->Id())
2445  {
2446  sout << "returning 1\n";
2447  }
2448  else {
2449  sout << "returning -1\n";
2450  }
2451  LOGPZ_DEBUG(loggerorient, sout.str())
2452  }
2453 #endif
2454 
2455  if(thisside.Element()->Id() < neighbour.Element()->Id())
2456  {
2457  return 1;
2458  }
2459  else {
2460  return -1;
2461  }
2462 
2463 }
2464 
2465 int ConjugateSide(TPZGeoEl *gel, int side, TPZStack<int> &allsides)
2466 {
2467  std::set<int> allside;
2468  allside.insert(&allsides[0],&allsides[0]+allsides.NElements());
2469  TPZStack<TPZGeoElSide> highsides;
2470  // the dimension of the conjugate side
2471  int targetdimension = gel->SideDimension(side)+1;
2472  // find all sides connected to side which have target dimension
2473  gel->AllHigherDimensionSides(side,targetdimension,highsides);
2474  int nhigh = highsides.NElements();
2475  int is;
2476  // find the side which is not contained in allsides
2477  for(is=0; is<nhigh; is++)
2478  {
2479  int highside = highsides[is].Side();
2480  if(allside.find(highside) == allside.end())
2481  {
2482  return highside;
2483  }
2484  }
2485  return -1;
2486 }
2487 
2488 void TPZGeoEl::GetHigherSubElements(TPZVec<TPZGeoEl*> &unrefinedSons)
2489 {
2490  int nsons = this->NSubElements();
2491  for(int s = 0; s < nsons; s++)
2492  {
2493  TPZGeoEl * son = this->SubElement(s);
2494  if(!son) continue;
2495 
2496  if(son->HasSubElement() == false)
2497  {
2498  int oldSize = unrefinedSons.NElements();
2499  unrefinedSons.Resize(oldSize+1, son);
2500  }
2501  else
2502  {
2503  son->GetHigherSubElements(unrefinedSons);
2504  }
2505  }
2506 }
2507 
2508 // Compute the permutation for an HDiv side
2509 void TPZGeoEl::HDivPermutation(int side, TPZVec<int> &permutegather)
2510 {
2511  int dimension = Dimension();
2512  int sidedimension = SideDimension(side);
2513 
2514  if(dimension != sidedimension+1)
2515  {
2516  std::stringstream sout;
2517  sout << "HDivPermutation called with wrong side parameter " << side;
2518 #ifdef LOG4CXX
2519  LOGPZ_ERROR(logger,sout.str())
2520 #endif
2521  cout << sout.str() << std::endl;
2522  }
2523  DebugStop();
2524 }
2525 
2526 void TPZGeoEl::GetNodeIndices( TPZVec<int64_t> &nodeindices ){
2527  const int nnodes = this->NNodes();
2528  nodeindices.Resize( nnodes );
2529  for(int i = 0; i < nnodes; i++){
2530  nodeindices[i] = this->NodeIndex(i);
2531  }
2532 }
2533 
2534 void TPZGeoEl::GetNodeIndices( std::set<int64_t> &nodeindices ){
2535  nodeindices.clear();
2536  const int nnodes = this->NNodes();
2537  for(int i = 0; i < nnodes; i++){
2538  nodeindices.insert(this->NodeIndex(i));
2539  }
2540 }
2541 
2542 
2543 void TPZGeoEl::NodesCoordinates(TPZFMatrix<REAL > &coordinates){
2544  coordinates.Resize(3,this->NNodes());
2545  coordinates.Zero();
2546  TPZManVector<REAL,3> co(3,0.0);
2547 
2548  for (int in = 0; in < this->NNodes(); in++) {
2549  TPZGeoNode inode = this->Node(in);
2550  inode.GetCoordinates(co);
2551  coordinates(0,in) = co[0];
2552  coordinates(1,in) = co[1];
2553  coordinates(2,in) = co[2];
2554  }
2555 
2556 }
2557 
2558 int TPZGeoEl::StaticClassId() {
2559  return Hash("TPZGeoEl");
2560 }
2561 
2562 int TPZGeoEl::ClassId() const {
2563  return StaticClassId();
2564 }
2565 
2566 TPZGeoNode* TPZGeoEl::NodePtr(int i) const {
2567  return &(fMesh->NodeVec()[NodeIndex(i)]);
2568 }
2569 
2570 TPZGeoNode& TPZGeoEl::Node(int i) const {
2571  return (fMesh->NodeVec()[NodeIndex(i)]);
2572 }
2573 
2574 TPZGeoEl* TPZGeoEl::Father() const {
2575  return (fFatherIndex == -1) ? 0 : Mesh()->ElementVec()[fFatherIndex];
2576 }
2577 
2578 TPZGeoNode* TPZGeoEl::SideNodePtr(int side, int nodenum) const {
2579  return &(fMesh->NodeVec()[SideNodeIndex(side, nodenum)]);
2580 }
TPZGeoEl::SideNodeIndex
virtual int64_t SideNodeIndex(int side, int nodenum) const =0
Returns the index of the nodenum node of side.
Normalize
void Normalize(TPZVec< REAL > &normlow, TPZVec< REAL > &normal)
Normalize normlow vector (in self)
Definition: pzgeoel.cpp:2093
TPZGeoEl::SmallerEdge
virtual REAL SmallerEdge()
Definition: pzgeoel.cpp:639
pzshapetetra.h
Contains TPZShapeTetra class which implements the shape functions of a tetrahedral element...
TPZVec::Print
void Print(std::ostream &out=std::cout)
Prints the structural information of the vector object to the output stream. This method will not p...
Definition: pzvec.h:482
fabs
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
TPZGeoEl::Father2
virtual TPZGeoElSide Father2(int side) const
Returns the father/side of the father which contains the side of the sub element. ...
Definition: pzgeoel.cpp:376
TPZGeoEl::Volume
REAL Volume()
Return the volume of the element.
Definition: pzgeoel.cpp:1052
TPZAdmChunkVector::AllocateNewElement
int AllocateNewElement()
Makes more room for new elements.
Definition: pzadmchunk.h:184
pzquad.h
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
TPZIntPoints::NPoints
virtual int NPoints() const =0
Returns number of points for the cubature rule related.
TPZGeoEl::NodePtr
TPZGeoNode * NodePtr(int i) const
Returns a pointer to the ith node of the element.
Definition: pzgeoel.cpp:2566
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
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.
ConjugateSide
int ConjugateSide(TPZGeoEl *gel, int side, TPZStack< int > &allsides)
Find a side which is not contained in allsides and whose dimension is one higher than the dimension o...
Definition: pzgeoel.cpp:2465
TPZGeoEl::QuadArea
static REAL QuadArea(TPZVec< TPZGeoNode *> &nodes)
Returns the area from a quadrilateral face.
Definition: pzgeoel.cpp:1032
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZFMatrix::GramSchmidt
void GramSchmidt(TPZFMatrix< TVar > &Orthog, TPZFMatrix< TVar > &TransfToOrthog)
This method implements a Gram Schimidt method. this = Orthog.TransfToOrthog.
Definition: pzfmatrix.cpp:341
TPZGeoEl::SetNeighbourInfo
virtual void SetNeighbourInfo(int side, TPZGeoElSide &neigh, TPZTransform< REAL > &trans)=0
TPZGeoEl::SetNeighbour
virtual void SetNeighbour(int side, const TPZGeoElSide &neighbour)=0
Fill in the data structure for the neighbouring information.
pzgnode.h
Contains declaration of TPZGeoNode class which defines a geometrical node.
TPZGeoElSide::NSubElements
int NSubElements()
return the number of element/side pairs which compose the current set of points
Definition: pzgeoelside.cpp:1010
TPZGeoNode::Coord
REAL Coord(int i) const
Returns i-th coordinate of the current node.
Definition: pzgnode.h:99
nsub
clarg::argInt nsub("-nsub", "number of substructs", 4)
dimension
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
TPZGeoEl::TriangleArea
static REAL TriangleArea(TPZVec< TPZGeoNode *> &nodes)
Returns the area from the triangular face.
Definition: pzgeoel.cpp:1009
TPZGeoEl::Distance
static REAL Distance(TPZVec< REAL > &centel, TPZVec< REAL > &centface)
Definition: pzgeoel.cpp:936
pzcompel.h
Contains declaration of TPZCompEl class which defines the interface of a computational element...
pzvec.h
Templated vector implementation.
TPZRefPattern::Element
TPZGeoEl * Element(int iel)
It returns the element number iel from the stack of elements of the geometric mesh.
Definition: TPZRefPattern.cpp:515
TPZGeoEl::IsGeoElMapped
virtual bool IsGeoElMapped() const
Returns if is a TPZGeoElMapped< T > element.
Definition: pzgeoel.h:282
TPZGeoElSide::Father2
TPZGeoElSide Father2() const
returns the father/side pair which contains the set of points associated with this/side ...
Definition: pzgeoelside.cpp:875
TPZGeoEl::BuildTransform2
virtual TPZTransform< REAL > BuildTransform2(int side, TPZGeoEl *father, TPZTransform< REAL > &t)
Returns the transformation which maps the parameter side of the element/side into the parameter spac...
Definition: pzgeoel.cpp:386
TPZGeoEl::Projection
TPZTransform< REAL > Projection(int side)
Compute the projection of the point within the interior of the element to the side of the element...
Definition: pzgeoel.cpp:1970
pzshape
groups all classes dedicated to the computation of shape functions
Definition: pzshapeextend.h:16
TPZGeoEl::fFatherIndex
int64_t fFatherIndex
Index of the element from which the element is a subelement.
Definition: pzgeoel.h:60
pzerror.h
Defines PZError.
TPZGeoEl::NCornerNodes
virtual int NCornerNodes() const =0
Returns the number of corner nodes of the element.
pzrefquad.h
Contains the TPZRefQuad class which implements the uniform refinement of a geometric quadrilateral el...
TPZGeoEl::Write
void Write(TPZStream &str, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzgeoel.cpp:1663
TPZGeoEl::SideNodePtr
virtual TPZGeoNode * SideNodePtr(int side, int nodenum) const
Returns the pointer to the nodenum node of side.
Definition: pzgeoel.cpp:2578
TPZGeoElSide::GetSubElements2
void GetSubElements2(TPZStack< TPZGeoElSide > &subelements)
build the list of element/side pairs which compose the current set of points
Definition: pzgeoelside.cpp:893
TPZGeoElSide
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
pzrefine::fatherside
static int fatherside[8][21]
Definition: pzrefprism.cpp:303
TPZGeoElSide::HasSubElement
int HasSubElement()
Return 1 if the element has subelements along side.
Definition: pzgeoelside.cpp:810
TPZChunkVector::NElements
int64_t NElements() const
Access method to query the number of elements of the vector.
Definition: TPZChunkVector.h:236
TPZGeoElSide::Neighbour
TPZGeoElSide Neighbour() const
Definition: pzgeoel.h:754
TPZGeoEl::Mesh
TPZGeoMesh * Mesh() const
Returns the mesh to which the element belongs.
Definition: pzgeoel.h:220
TPZGeoEl::NSides
virtual int NSides() const =0
Returns the number of connectivities of the element.
TPZGeoEl::Divide
virtual void Divide(TPZVec< TPZGeoEl *> &pv)
Divides the element and puts the resulting elements in the vector.
Definition: pzgeoel.h:742
TPZPersistenceManager::GetInstance
static TPZSavable * GetInstance(const int64_t &objId)
Definition: TPZPersistenceManager.cpp:397
TPZGeoEl::NodeIndex
virtual int64_t NodeIndex(int i) const =0
Returns the index of the ith node the index is the location of the node in the nodevector of the mesh...
TPZGeoElSide::CenterPoint
void CenterPoint(TPZVec< REAL > &center) const
return the coordinates of the center in master element space (associated with the side) ...
Definition: pzgeoelside.cpp:414
TPZGeoEl::SideDimension
virtual int SideDimension(int side) const =0
Return the dimension of side.
TPZGeoElSide::SideTransform3
void SideTransform3(TPZGeoElSide neighbour, TPZTransform<> &t)
Accumulates the transformations from the current element/side to the neighbour/side.
Definition: pzgeoelside.cpp:601
TPZGeoEl::fMatId
int fMatId
Material index.
Definition: pzgeoel.h:54
TPZGeoEl::WhichSide
int WhichSide(TPZVec< int64_t > &SideNodeIds)
Returns the side number which is connected to the SideNodes returns -1 if no side is found...
Definition: pzgeoel.cpp:165
TPZGeoEl::CheckSubelDataStructure
void CheckSubelDataStructure()
Checks the structure of the Father() and GetSubelement2()
Definition: pzgeoel.cpp:503
pzshapelinear.h
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
pzrefpyram.h
Contains the TPZRefPyramid class which implements the uniform refinement of a geometric hexahedral el...
TPZGeoMesh::SetElementIdUsed
void SetElementIdUsed(int64_t id)
Indicates that an element with id was created.
Definition: pzgmesh.h:111
TPZGeoElSide::BuildConnectivities
static void BuildConnectivities(TPZVec< TPZGeoElSide > &elvec, TPZVec< TPZGeoElSide > &neighvec)
Definition: pzgeoelside.cpp:1017
TPZIntPoints
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
error
void error(char *string)
Definition: testShape.cc:7
TPZGeoEl::ComputeParamTrans
TPZTransform< REAL > ComputeParamTrans(TPZGeoEl *fat, int fatside, int sideson)
Definition: pzgeoel.cpp:845
TPZGeoEl::CenterPoint
virtual void CenterPoint(int side, TPZVec< REAL > &masscent) const =0
It returns the coordinates from the center of the side of the element in the element coordinate space...
TPZGeoEl::Read
void Read(TPZStream &str, void *context) override
read objects from the stream
Definition: pzgeoel.cpp:1652
TPZGeoEl::EldestAncestor
TPZGeoEl * EldestAncestor() const
Definition: pzgeoel.cpp:1981
pzshapecube.h
Contains TPZShapeCube class which implements the shape functions of a hexaedral element.
NormalVector
void NormalVector(TPZGeoElSide &LC, TPZGeoElSide &LS, TPZVec< REAL > &normal)
Computes the normal vector goes from the center of the conjugate side (LC) to the center of the LS si...
Definition: pzgeoel.cpp:1999
TPZGeoEl::SubElement
virtual TPZGeoEl * SubElement(int is) const =0
Returns a pointer to the subelement is.
TPZGeoEl::CreateSideIntegrationRule
virtual TPZIntPoints * CreateSideIntegrationRule(int side, int order)=0
Creates an integration rule for the topology of the corresponding side and able to integrate a polyno...
pztopology::highsides
static int highsides[27][7]
For each side was stored the sides connected with it but of the higher dimension. ...
Definition: tpzcube.cpp:110
TPZCompMesh::CreateCompEl
TPZCompEl * CreateCompEl(TPZGeoEl *gel, int64_t &index)
Create a computational element based on the geometric element.
Definition: pzcmesh.h:462
TPZGeoEl::fId
int64_t fId
Traditional element number or element id.
Definition: pzgeoel.h:51
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZVec::Resize
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
pzrefprism.h
Contains the TPZRefPrism class which implements the uniform refinement of a geometric prism element...
TPZGeoEl::operator<<
friend std::ostream & operator<<(std::ostream &s, TPZGeoEl &el)
Equivalent to Print.
Definition: pzgeoel.cpp:314
TPZGeoEl::GetAllSiblings
virtual void GetAllSiblings(TPZStack< TPZGeoEl *> &unrefinedSons)
[deprecated] use YoungestChildren
Definition: pzgeoel.cpp:410
TPZGeoEl::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzgeoel.cpp:2562
TPZGeoLinear.h
TPZGeoMesh::CreateUniqueElementId
int64_t CreateUniqueElementId()
Returns ++fElementMaxId.
Definition: pzgmesh.h:117
pzgmesh.h
Contains declaration of TPZMesh class which defines a geometrical mesh and contains a corresponding l...
TPZGeoEl::TransformSonToFather
void TransformSonToFather(TPZGeoEl *ancestor, TPZVec< REAL > &ksiSon, TPZVec< REAL > &ksiAncestor)
Compute the map of a paramenter point in the subelement to a parameter point in the super element...
Definition: pzgeoel.cpp:816
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
TPZGeoEl::Level
int Level()
Returns the number of ancestors.
Definition: pzgeoel.cpp:319
TPZGeoEl::CreateBCCompEl
virtual TPZCompEl * CreateBCCompEl(int side, int bc, TPZCompMesh &cmesh)
Method which creates a computational boundary condition element based on the current geometric elemen...
Definition: pzgeoel.cpp:1092
TPZGeoEl::ElementExists
int ElementExists(TPZGeoEl *elem, int64_t id)
To test continuity.
Definition: pzgeoel.cpp:361
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZGeoMesh::Print
virtual void Print(std::ostream &out=std::cout) const
Print the information of the grid to an ostream.
Definition: pzgmesh.cpp:130
TPZGeoEl::WhichSubel
int WhichSubel() const
Returns the son number of the sub element gel.
Definition: pzgeoel.cpp:448
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
Max
const T & Max(const T &a, const T &b)
Returns the maximum value between a and b.
Definition: pzreal.h:724
TPZGeoEl::Index
int64_t Index() const
Returns the index of the element within the element vector of the mesh.
Definition: pzgeoel.h:730
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
pzgeotetrahedra.h
Contains the TPZGeoTetrahedra class which implements the geometry of a tetrahedral element...
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZGeoEl::FatherSide
virtual int FatherSide(int side, int son)
Definition: pzgeoel.cpp:381
TPZIntPoints::SetOrder
virtual void SetOrder(TPZVec< int > &ord, int type=0)=0
Sets the order of the cubature rule.
TPZGeoEl::MidSideNodeIndices
virtual void MidSideNodeIndices(int side, TPZVec< int64_t > &indices) const
Returns the midside node indices along a side of the element.
Definition: pzgeoel.cpp:1137
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZGeoElSide::X
void X(TPZVec< REAL > &loc, TPZVec< REAL > &result) const
X coordinate of a point loc of the side.
Definition: pzgeoelside.cpp:117
TPZGeoEl::ShapePhi1d
void ShapePhi1d(REAL x, int num, TPZFMatrix< REAL > &phi)
Computes the values of the "num" one dimensional shapefunctions at point x, using lagrangian interpol...
Definition: pzgeoel.cpp:149
TPZGeoEl::GetHigherSubElements
virtual void GetHigherSubElements(TPZVec< TPZGeoEl *> &unrefinedSons)
Definition: pzgeoel.cpp:2488
ELU
Definition: pzmatrix.h:52
TPZGeoEl::fReference
TPZCompEl * fReference
Reference to the element currently loaded. Pointer is given as this->Mesh()->Reference()->ElementVec(...
Definition: pzgeoel.h:57
TPZGeoEl::fNumInterfaces
int fNumInterfaces
A counter to indicate how many interface elements are pointing to it.
Definition: pzgeoel.h:69
TPZGeoElSide::RemoveConnectivity
void RemoveConnectivity()
Remove the element from the connectivity loop.
Definition: pzgeoelside.cpp:313
TPZGeoEl::GetSubElements2
virtual void GetSubElements2(int side, TPZStack< TPZGeoElSide > &subel) const
This method will return a partition of the side of the current element as the union of sub elements/...
Definition: pzgeoel.cpp:392
TPZGeoEl::Jacobian
void Jacobian(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &jac, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Compute a decomposition of the gradient of the mapping function, as a rotation matrix (Jacobian) and ...
Definition: pzgeoel.cpp:1144
TPZGeoMesh::NodeVec
TPZAdmChunkVector< TPZGeoNode > & NodeVec()
Definition: pzgmesh.h:140
TPZGeoEl::Reference
TPZCompEl * Reference() const
Return a pointer to the element referenced by the geometric element.
Definition: pzgeoel.h:750
TPZGeoEl::BuildBlendConnectivity
void BuildBlendConnectivity()
Set connectivity information elements with blend geometric map.
Definition: pzgeoel.cpp:1959
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZGeoEl::NormalOrientation
int NormalOrientation(int side)
Determine the orientation of the normal vector comparing the ids of the neighbouring elements...
Definition: pzgeoel.cpp:2370
TPZGeoEl::fIndex
int64_t fIndex
Index of the element in the element vector.
Definition: pzgeoel.h:63
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZGeoEl::CreateBCGeoEl
virtual TPZGeoEl * CreateBCGeoEl(int side, int bc)=0
Method which creates a geometric element on the side of an existing element.
TPZGeoEl::InitializeNeighbours
void InitializeNeighbours()
To be used after the buid connectivity. If some neighbour isn&#39;t initialized.
Definition: pzgeoel.cpp:1117
tpzgeoelrefpattern.h
Contains declaration of TPZGeoElRefPattern class which implements a generic geometric element which i...
pzrefine
Groups all classes which model the h refinement These classes are used as template arguments of...
Definition: pzrefpoint.cpp:15
co
REAL co[8][3]
Coordinates of the eight nodes.
Definition: tpzgensubstruct.cpp:56
TPZCompEl::Index
int64_t Index() const
Returns element index of the mesh fELementVec list.
Definition: pzcompel.h:821
Min
const T & Min(const T &a, const T &b)
Returns the minimum value between a and b.
Definition: pzreal.h:729
pzshapepoint.h
Contains TPZShapePoint class which implements the shape function associated with a point...
TPZGeoEl::Neighbour
virtual TPZGeoElSide Neighbour(int side)=0
Returns a pointer to the neighbour and the neighbourside along side of the current element...
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZGeoEl::JacobianXYZ
void JacobianXYZ(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &jac, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Compute a decomposition of the gradient of the mapping function, as a rotation matrix (Jacobian) and ...
Definition: pzgeoel.cpp:1159
pzcmesh.h
Contains declaration of TPZCompMesh class which is a repository for computational elements...
TPZRefPattern.h
Contains the TPZRefPattern class which defines the topology of the current refinement pattern to a me...
TPZGeoEl::Node
TPZGeoNode & Node(int i) const
Returns the ith node of the element.
Definition: pzgeoel.cpp:2570
TPZGeoCube.h
Contains the TPZGeoCube class which implements the geometry of hexahedra element. ...
TPZGeoEl::GetNodeIndices
void GetNodeIndices(TPZVec< int64_t > &nodeindices)
Definition: pzgeoel.cpp:2526
TPZGeoEl::NSideNodes
virtual int NSideNodes(int side) const =0
Returns the number of nodes for a particular side.
pzstack.h
A simple stack.
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
pzgeopoint.h
Contains the TPZGeoPoint class which implements the geometry of a point element or 0-D element...
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZGeoEl::HasSubElement
virtual int HasSubElement() const =0
Return 1 if the element has subelements.
TPZTransform::Multiply
TPZTransform< T > Multiply(TPZTransform< T > &right)
Multiply the transformation object (to the right) with right (Multiplying matrices) ...
Definition: pztrnsform.cpp:112
TPZGeoEl::IsInParametricDomain
virtual bool IsInParametricDomain(const TPZVec< REAL > &pt, REAL tol=1.e-2)=0
Verifies if the parametric point pt is in the element parametric domain.
TPZGeoEl::NNodes
virtual int NNodes() const =0
Returns the number of nodes of the element.
pzshapeprism.h
Contains TPZShapePrism class which implements the shape functions of a prism element.
TPZGeoEl::IsLinearMapping
virtual bool IsLinearMapping() const
Definition: pzgeoel.h:268
TPZGeoElSide::Element
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
TPZGeoEl::NeighbourExists
int NeighbourExists(int side, const TPZGeoElSide &gel)
Returns 1 if gel is a neighbour of the element along side.
Definition: pzgeoel.cpp:226
TPZGeoEl::PrintTopologicalInfo
virtual void PrintTopologicalInfo(std::ostream &out=std::cout)
Prints the coordinates of all nodes (geometric)
Definition: pzgeoel.cpp:300
TPZGeoEl::HDivPermutation
virtual void HDivPermutation(int side, TPZVec< int > &permutegather)
Computes the permutation for an HDiv side.
Definition: pzgeoel.cpp:2509
pzreftriangle.h
Contains the TPZRefTriangle class which implements the uniform refinement of a geometric triangular e...
TPZGeoEl::GradX
virtual void GradX(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &gradx) const =0
Return the gradient of the transformation at the given coordinate.
TPZRefPattern::fId
int fId
id of the refinement pattern
Definition: TPZRefPattern.h:124
TPZGeoNode::Print
void Print(std::ostream &out=std::cout)
Print the node data into out.
Definition: pzgnode.cpp:88
TPZGeoEl::~TPZGeoEl
virtual ~TPZGeoEl()
Destructor.
Definition: pzgeoel.cpp:49
pzreftetrahedra.h
Contains the TPZRefTetrahedra class which implements the uniform refinement of a geometric tetrahedra...
TPZGeoEl::Dimension
virtual int Dimension() const =0
Returns the dimension of the element.
TPZGeoEl::SideToSideTransform
virtual TPZTransform< REAL > SideToSideTransform(int sidefrom, int sideto)=0
Compute the transformation between the master element space of one side of an element to the master e...
TPZRefPattern
Defines the topology of the current refinement pattern to a mesh. Refine.
Definition: TPZRefPattern.h:77
dist
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
TPZGeoEl::AllHigherDimensionSides
virtual void AllHigherDimensionSides(int side, int targetdimension, TPZStack< TPZGeoElSide > &elsides)=0
TPZGeoElSide::IsRelative
bool IsRelative(TPZGeoElSide other)
Checks whether other is a relative (son or ancestor) of this.
Definition: pzgeoelside.cpp:111
TPZGeoEl::Id
int64_t Id() const
Returns the Id of the element.
Definition: pzgeoel.h:223
TPZFMatrix::Read
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: pzfmatrix.cpp:2205
TPZGeoEl::ComputeXInverse
bool ComputeXInverse(TPZVec< REAL > &XD, TPZVec< REAL > &ksi, REAL Tol)
Computes the XInverse and returns true if ksi belongs to master element domain.
Definition: pzgeoel.cpp:686
TPZGeoNode
Implements a geometric node in the pz environment. Geometry.
Definition: pzgnode.h:31
TPZGeoEl::X
virtual void X(TPZVec< REAL > &qsi, TPZVec< REAL > &result) const =0
Return the coordinate in real space of the point coordinate in the master element space...
TPZGeoEl::RefElVolume
virtual REAL RefElVolume()=0
Volume of the master element.
TPZGeoEl::SideArea
virtual REAL SideArea(int side)
Returns the area from the face.
Definition: pzgeoel.cpp:1064
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
TPZGeoEl::IsGeoBlendEl
virtual bool IsGeoBlendEl() const
Definition: pzgeoel.h:275
TPZStack
This class implements a stack object. Utility.
Definition: pzcheckmesh.h:14
TPZGeoElSide::Area
REAL Area()
Area associated with the side.
Definition: pzgeoelside.cpp:259
TPZGeoElSide::Dimension
int Dimension() const
the dimension associated with the element/side
Definition: pzgeoelside.cpp:593
TPZGeoElSide::ResetBlendConnectivity
bool ResetBlendConnectivity(const int64_t &index)
Definition: pzgeoelside.cpp:338
TPZGeoElSide::Jacobian
void Jacobian(TPZVec< REAL > &param, TPZFMatrix< REAL > &jacobian, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Jacobian associated with the side of the element.
Definition: pzgeoelside.cpp:189
TPZCompMesh
Implements computational mesh. Computational Mesh.
Definition: pzcmesh.h:47
TPZGeoEl::CharacteristicSize
virtual REAL CharacteristicSize()
Computes the set of normals for defining HDiv approximation spaces.
Definition: pzgeoel.cpp:604
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
pztopology::sidedimension
static int sidedimension[27]
Vector of the dimension for each side.
Definition: tpzcube.cpp:100
TPZGeoEl::RemoveConnectivities
void RemoveConnectivities()
It removes the connectivities of the element.
Definition: pzgeoel.cpp:1098
TPZGeoEl::TPZGeoEl
TPZGeoEl()
Definition: pzgeoel.h:173
TPZGeoEl::GetTransformId2dT
int GetTransformId2dT(TPZVec< int > &idfrom, TPZVec< int > &idto)
Get the transform id the face to face.
Definition: pzgeoel.cpp:347
TPZGeoElSide::Side
int Side() const
Definition: pzgeoelside.h:169
TPZRefLinear.h
Contains the TPZRefLinear class which implements the uniform refinement of a geometric linear element...
TPZGeoEl::GetRefPattern
virtual TPZAutoPointer< TPZRefPattern > GetRefPattern() const
Returns the refinement pattern associated with the element.
Definition: pzgeoel.cpp:1716
TPZGeoEl::gGlobalAxes
static TPZFMatrix< REAL > gGlobalAxes
3x3 unit matrix to be copied to the axes if the geometric element does not have a particular orientat...
Definition: pzgeoel.h:66
TPZGeoEl::Shape1d
void Shape1d(REAL x, int num, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
Computes the values of the "num" one dimensional shapefunctions and derivatives at point x...
Definition: pzgeoel.cpp:128
TPZGeoEl::SetRefPattern
virtual void SetRefPattern(TPZAutoPointer< TPZRefPattern >)
Defines the refinement pattern. It&#39;s used only in TPZGeoElRefPattern objects.
Definition: pzgeoel.cpp:1647
TPZGeoEl::NodesCoordinates
void NodesCoordinates(TPZFMatrix< REAL > &cooridnates)
Computes nodes coordinates of the element.
Definition: pzgeoel.cpp:2543
TPZGeoElSide::Exists
int Exists() const
Definition: pzgeoelside.h:175
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZGeoEl::fMesh
TPZGeoMesh * fMesh
Pointer to the mesh to which the element belongs.
Definition: pzgeoel.h:48
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
pzshapepiram.h
Contains TPZShapePiram class which implements the shape functions of a pyramid element.
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
pzgeopyramid.h
Contains the TPZGeoPyramid class which implements the geometry of pyramid element.
pztrnsform.h
Contains the TPZTransform<> class which implements an affine transformation between points in paramet...
TPZGeoEl::VerifyNodeCoordinates
bool VerifyNodeCoordinates(REAL tol=1e-1)
Verify coordinate of element nodes checking if they are coincident to the X mapping of the corner nod...
Definition: pzgeoel.cpp:1722
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
rdt.values
def values
Definition: rdt.py:119
TPZTransform::SetMatrix
void SetMatrix(TPZFMatrix< T > &mult, TPZFMatrix< T > &sum)
Sets the transformation matrices.
Definition: pztrnsform.cpp:96
pzgeoprism.h
Contains the TPZGeoPrism class which implements the geometry of a prism element.
TPZGeoEl::ElementRadius
virtual REAL ElementRadius()
Definition: pzgeoel.cpp:949
TPZGeoNode::GetCoordinates
void GetCoordinates(TPZVec< REAL > &co)
Fill the coordinates of the node.
Definition: pzgnode.cpp:83
TPZGeoEl::Print
virtual void Print(std::ostream &out=std::cout)
Print all relevant data of the element to cout.
Definition: pzgeoel.cpp:238
pzgeom
Groups all classes which model the geometry.
Definition: pzgeopoint.cpp:18
TPZGeoEl::StaticClassId
static int StaticClassId()
Definition: pzgeoel.cpp:2558
TPZFMatrix::Write
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzfmatrix.cpp:2225
ZeroTolerance
REAL ZeroTolerance()
Returns the tolerance to Zero value. Actually: .
Definition: pzreal.h:633
TPZGeoEl::SetNeighbourForBlending
void SetNeighbourForBlending(int side)
TPZGeoBlend need to find a non-linear neighbour.
Definition: pzgeoel.cpp:1913
TPZSavable
This class defines the interface to save and restore objects from TPZStream objects. Persistency.
Definition: TPZSavable.h:67
TPZCompEl
Defines the interface of a computational element. Computational Element.
Definition: pzcompel.h:59
TPZGeoEl::SetSubElementConnectivities
virtual void SetSubElementConnectivities()
Initializes the external connectivities of the subelements.
Definition: pzgeoel.cpp:563
TPZGeoEl::NSubElements
virtual int NSubElements() const =0
Returns the number of subelements of the element independent of the fact whether the element has alr...
TPZGeoEl::Father
TPZGeoEl * Father() const
Computes the normal vectors needed for forming HDiv vector valued shape functions.
Definition: pzgeoel.cpp:2574
pzshapequad.h
Contains TPZShapeQuad class which implements the shape functions of a quadrilateral element...
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
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
TPZGeoEl::GetTransformId2dQ
int GetTransformId2dQ(TPZVec< int > &idfrom, TPZVec< int > &idto)
Get the transform id the face to face.
Definition: pzgeoel.cpp:330
TPZGeoElSide::SideNodeIndex
int64_t SideNodeIndex(int nodenum) const
Returns the index of the nodenum node of side.
Definition: pzgeoelside.cpp:821
TPZRefCube.h
Contains the TPZRefCube class which implements the uniform refinement of a geometric hexahedral eleme...
pzrefpoint.h
Contains the TPZRefPoint class which implements the uniform refinement of a geometric point element...
TPZPersistenceManager::WritePointer
static void WritePointer(const TPZSavable *obj, TPZStream *stream)
Definition: TPZPersistenceManager.cpp:182
pzshapetriang.h
Contains TPZShapeTriang class which implements the shape functions of a triangular element...
TPZTransform::Apply
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
TPZIntPoints::Point
virtual void Point(int i, TPZVec< REAL > &pos, REAL &w) const =0
Returns i-th point at master element and related weight.
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZGeoMesh::ElementVec
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138
TPZGeoEl::MidSideNodeIndex
virtual void MidSideNodeIndex(int side, int64_t &index) const =0
Returns the midside node index along a side of the element.
TPZGeoEl::YoungestChildren
virtual void YoungestChildren(TPZStack< TPZGeoEl *> &unrefinedSons)
This method will return all children at the bottom of the refinement tree of the element. i.e. all children that have no subelements.
Definition: pzgeoel.cpp:431
Generated on Sun Aug 16 2020 11:17:06 for NeoPZ by   doxygen 1.8.13