NeoPZ
pzgeoelside.cpp
Go to the documentation of this file.
1 
6 #include "pzgeoelside.h"
7 #include "pzgeoel.h"
8 #include "pztrnsform.h"
9 #include "pzstack.h"
10 #include "pzvec_extras.h"
11 #include "pzquad.h"
12 #include "pzshapequad.h"
13 #include "pzshapetriang.h"
14 #include "pzcompel.h"
15 #include "pzintel.h"
16 #include "pznumeric.h"
17 
18 #include "pzmultiphysicscompel.h"
19 
20 using namespace pzshape;
21 using namespace std;
22 
23 #include "pzlog.h"
24 
25 #include <sstream>
26 
27 #ifdef LOG4CXX
28 static LoggerPtr logger(Logger::getLogger("pz.mesh.tpzgeoelside"));
29 #endif
30 
31 // Implementation of the TPZGeoElSideIndex methods
32 
34  if (gel) this->fGeoElIndex = gel->Index();
35  else this->fGeoElIndex = -1;
36  this->fSide = side;
37 }
38 
40  TPZGeoEl * gel = side.Element();
41  if (gel) this->fGeoElIndex = gel->Index();
42  else this->fGeoElIndex = -1;
43  this->fSide = side.Side();
44 }
45 
47  if (geoel) this->fGeoElIndex = geoel->Index();
48  else this->fGeoElIndex = -1;
49 }
50 
52  return Hash("TPZGeoElSideIndex");
53 }
54 
55 void TPZGeoElSideIndex::Read(TPZStream& buf, void* context) { //ok
56  buf.Read(&fGeoElIndex);
57  buf.Read(&fSide);
58 }
59 
60 void TPZGeoElSideIndex::Write(TPZStream& buf, int withclassid) const { //ok
61  buf.Write(&fGeoElIndex);
62  buf.Write(&fSide);
63 }
64 
65 
66 // Implementation of the TPZGeoElSide methods
67 
68 TPZGeoElSide::TPZGeoElSide(TPZGeoEl *gel, std::set<int64_t> &sideCornerNodes)
69 {
70  fGeoEl = 0; fSide = -1;
71 
72  std::set<int64_t> nodes;
73  for(int s = 0; s < gel->NSides(); s++)
74  {
75  nodes.clear();
76  TPZGeoElSide actualSide(gel, s);
77  if(actualSide.NSideNodes() != (int)sideCornerNodes.size())
78  {
79  continue;
80  }
81  for(int as = 0; as < actualSide.NSideNodes(); as++)
82  {
83  nodes.insert(actualSide.SideNodeIndex(as));
84  }
85  if(nodes == sideCornerNodes)
86  {
87  fGeoEl = gel;
88  fSide = s;
89 
90  return;
91  }
92  }
93 }
94 
96  if(*this == other) return true;
97  TPZGeoElSide father = this->Father2();
98  if(father.Element()){
99  if(father.Element() == other.Element()){
100  return true;
101  }
102  else{
103  if(father.IsAncestor(other)){
104  return true;
105  }
106  }
107  }
108  return false;
109 }
110 
112  if( this->IsAncestor(other) ) return true;
113  if( other.IsAncestor(*this) ) return true;
114  return false;
115 }
116 
117 void TPZGeoElSide::X(TPZVec< REAL > &loc, TPZVec< REAL > &result) const {
118 
119  TPZManVector< REAL,3 > locElement(fGeoEl->Dimension(), 0.);
120  result.Resize(3);
121  QsiElement(loc, locElement);
122  fGeoEl->X(locElement, result);
123 }
124 
126 void TPZGeoElSide::QsiElement(TPZVec< REAL > &qsi_side, TPZVec< REAL > &qsi_element) const{
127 
128  qsi_element.Resize(fGeoEl->Dimension());
129  TPZTransform<> ElementDim = fGeoEl->SideToSideTransform(fSide, fGeoEl->NSides()-1);
130  ElementDim.Apply(qsi_side, qsi_element);
131 }
132 
135 
136 #ifdef PZDEBUG
137  if(!fGeoEl) return;
138 #endif
139 
140  int dim = fGeoEl->Dimension();
141  TPZFNMatrix<9,REAL> gradx_vol(3,dim);
142 
143  TPZManVector< REAL,3 > locElement(dim, 0.);
144 
145  TPZTransform<> Transformation = fGeoEl->SideToSideTransform(fSide, fGeoEl->NSides()-1);
146 
147  Transformation.Apply(loc, locElement);
148  TPZFNMatrix<9,REAL> trans_mult(Transformation.Mult().Cols(),Transformation.Mult().Rows());
149  Transformation.Mult().Transpose(&trans_mult);
150  fGeoEl->GradX(locElement, gradx_vol);
151  gradx_vol.Multiply(Transformation.Mult(), gradx);
152 
153 }
154 
155 #ifdef _AUTODIFF
156 
157 void TPZGeoElSide::X(TPZVec< Fad<REAL> > &loc, TPZVec< Fad<REAL> > &result) const
158 {
159  TPZManVector<Fad<REAL>,3 > locElement(fGeoEl->Dimension(), 0.);
160  result.Resize(3);
161 
162  TPZTransform<> ElementDimR = fGeoEl->SideToSideTransform(fSide, fGeoEl->NSides()-1);
163  TPZTransform<Fad<REAL> > ElementDim;
164  ElementDim.CopyFrom(ElementDimR);
165 
166  ElementDim.Apply(loc, locElement);
167 
168  fGeoEl->X(locElement, result);
169 
170 }
171 
173 void TPZGeoElSide::GradX(TPZVec< Fad<REAL> > &loc, TPZFMatrix< Fad<REAL> > &gradx) const{
174 
175  TPZManVector< Fad<REAL> ,3 > locElement(fGeoEl->Dimension(), 0.);
176  gradx.Resize(3,fGeoEl->Dimension());
177 
178  TPZTransform<> ElementDimR = fGeoEl->SideToSideTransform(fSide, fGeoEl->NSides()-1);
179  TPZTransform<Fad<REAL> > ElementDim;
180  ElementDim.CopyFrom(ElementDimR);
181  ElementDim.Apply(loc, locElement);
182  fGeoEl->GradX(locElement, gradx);
183 }
184 
185 #endif
186 
187 
188 
189 void TPZGeoElSide::Jacobian(TPZVec<REAL> &param,TPZFMatrix<REAL> &jacobian,TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv) const {
190 
191  if(!fGeoEl) return;
192  int DIM = fGeoEl->Dimension();
193 
194  TPZTransform<> ThisTransf = fGeoEl->SideToSideTransform(fSide, fGeoEl->NSides()-1);
195  TPZManVector< REAL,3 > paramElement(DIM,0.);
196 
197  ThisTransf.Apply(param,paramElement);
198  REAL DetElement;
199 
200  TPZFNMatrix<9> JCn(3,DIM,0.), JacElement(DIM,DIM,0.), AxesElement(DIM,3,0.), Temp(3,DIM,0.), InvElement(DIM,DIM,0.), axest(3,DIM);
201 
202  fGeoEl->Jacobian(paramElement,JacElement,AxesElement,DetElement,InvElement);
203 
204  AxesElement.Transpose();
205  AxesElement.Multiply(JacElement,JCn);
206  JCn.Multiply(ThisTransf.Mult(),Temp);
207 
208  Temp.GramSchmidt(axest,jacobian);
209  axest.Transpose(&axes);
210  jacinv.Resize(jacobian.Rows(),jacobian.Cols());
211  if(axes.Rows() == 1)
212  {
213  detjac = jacobian(0,0);
214  if(detjac)
215  {
216  jacinv(0,0) = 1./jacobian(0,0);
217  }
218  else
219  {
220  jacinv(0,0) = 0.;
221  }
222  }
223  if(axes.Rows() == 2) {
224  detjac = jacobian(0,0)*jacobian(1,1) - jacobian(1,0)*jacobian(0,1);
225  if(detjac)
226  {
227  jacinv(0,0) = jacobian(1,1)/detjac; jacinv(0,1) = -jacobian(0,1)/detjac;
228  jacinv(1,0) = -jacobian(1,0)/detjac; jacinv(1,1) = jacobian(0,0)/detjac;
229  }
230  else
231  {
232  jacinv.Zero();
233  }
234  }
235  if(axes.Rows() == 3) {
236  detjac = -(jacobian(0,2)*jacobian(1,1)*jacobian(2,0)) + jacobian(0,1)*jacobian(1,2)*jacobian(2,0) + jacobian(0,2)*jacobian(1,0)*jacobian(2,1) - jacobian(0,0)*jacobian(1,2)*jacobian(2,1) - jacobian(0,1)*jacobian(1,0)*jacobian(2,2) + jacobian(0,0)*jacobian(1,1)*jacobian(2,2);
237  jacinv(0,0) = (-(jacobian(1,2)*jacobian(2,1)) + jacobian(1,1)*jacobian(2,2))/detjac;
238  jacinv(0,1) = (jacobian(0,2)*jacobian(2,1) - jacobian(0,1)*jacobian(2,2))/detjac;
239  jacinv(0,2) = (-(jacobian(0,2)*jacobian(1,1)) + jacobian(0,1)*jacobian(1,2))/detjac;
240  jacinv(1,0) = (jacobian(1,2)*jacobian(2,0) - jacobian(1,0)*jacobian(2,2))/detjac;
241  jacinv(1,1) = (-(jacobian(0,2)*jacobian(2,0)) + jacobian(0,0)*jacobian(2,2))/detjac;
242  jacinv(1,2) = (jacobian(0,2)*jacobian(1,0) - jacobian(0,0)*jacobian(1,2))/detjac;
243  jacinv(2,0) = (-(jacobian(1,1)*jacobian(2,0)) + jacobian(1,0)*jacobian(2,1))/detjac;
244  jacinv(2,1) = (jacobian(0,1)*jacobian(2,0) - jacobian(0,0)*jacobian(2,1))/detjac;
245  jacinv(2,2) = (-(jacobian(0,1)*jacobian(1,0)) + jacobian(0,0)*jacobian(1,1))/detjac;
246  }
247 }
248 
251 {
252  TPZStack<int> lower;
253  fGeoEl->LowerDimensionSides(fSide,lower);
254  return lower.NElements()+1;
255 }
256 
257 
260 {
261  TPZManVector<REAL,3> sideparam(Dimension(),0);
262  REAL detjac;
263  TPZFNMatrix<9> jacinv(3,3),jacobian(3,3),axes(3,3);
264  //supondo jacobiano constante: X linear
265  CenterPoint(sideparam);
266  Jacobian(sideparam,jacobian,axes,detjac,jacinv);
267  TPZIntPoints *intrule = fGeoEl->CreateSideIntegrationRule(fSide, 0);
268  REAL RefElVolume = 0.;
269  int np = intrule->NPoints();
270  TPZManVector<REAL,3> points(Dimension());
271  REAL weight;
272  for (int ip=0; ip<np; ip++) {
273  intrule->Point(ip, points, weight);
274  RefElVolume += weight;
275  }
276  return (RefElVolume*detjac);//RefElVolume(): volume do elemento de refer�ncia
277 
278 }
279 
280 
282 {
283  int nneighbours = 0;
284  TPZGeoElSide neigh = this->Neighbour();
285 
286  while(neigh.Element() != this->Element())
287  {
288  nneighbours++;
289  neigh = neigh.Neighbour();
290  }
291 
292  return nneighbours;
293 }
294 
295 
297 {
298  int nneighbours = 0;
299  TPZGeoElSide neigh = this->Neighbour();
300 
301  while(neigh.Element() != this->Element())
302  {
303  if(neigh.Element() != thisElem)
304  {
305  nneighbours++;
306  }
307  neigh = neigh.Neighbour();
308  }
309 
310  return nneighbours;
311 }
312 
314 
315  if(!Exists()) return;
316  if(fSide < 0 || fSide >= fGeoEl->NSides()) {
317  PZError << "TPZGeoElSide::SetConnectivity Index out of bound\n";
318  }
319  //it removes the connectivity of the cycle where this inserted one:
320  //neighpre->this->neighpos => neighpre->neighpos
321  TPZGeoElSide neighpre,neigh = Neighbour();
322  if(neigh.Element() == NULL || neigh.Side() == -1){
323  PZError << "TPZGeoElSide::SetConnectivity trying to remove null or inexistent connection";
324  }
325  TPZGeoElSide neighpos = neigh;
326  while(neigh.Exists() && neigh != *this){
327  neighpre = neigh;
328  neigh = neigh.Neighbour();
329  }
330  if(neigh == *this){
331  this->SetNeighbour(TPZGeoElSide());
332  if (neighpre.Exists()) neighpre.SetNeighbour(neighpos);
333  } else {
334  PZError << "TPZGeoElSide::SetConnectivity neighbourhood cycle error";
335  }
336 }
337 
338 bool TPZGeoElSide::ResetBlendConnectivity(const int64_t &index){
339  return fGeoEl->ResetBlendConnectivity(fSide, index);
340 }
341 using namespace std;
342 
347 {
348  if(!fGeoEl || !neighbour.fGeoEl) return;
349 
350  TPZStack<int> mylowerdimension, neighbourlowerdimension;
351  fGeoEl->LowerDimensionSides(fSide,mylowerdimension);
352  neighbour.fGeoEl->LowerDimensionSides(neighbour.fSide,neighbourlowerdimension);
353  SetConnectivity(neighbour);
354  int ns = mylowerdimension.NElements();
355  int is;
356  for(is=0; is<ns; is++)
357  {
358  TPZGeoElSide myl(fGeoEl,mylowerdimension[is]);
359  TPZGeoElSide neighl(neighbour.fGeoEl,neighbourlowerdimension[is]);
360  myl.SetConnectivity(neighl);
361  }
362 
363 }
364 
365 void TPZGeoElSide::SetConnectivity(const TPZGeoElSide &neighbour) const{
366 
367  if(!Exists()) return;
368  if(fSide >= fGeoEl->NSides()) {
369  PZError << "ERROR(TPZGeoElSide::SetConnectivity)-> Index greater than number of sides.\n";
370  PZError << " fNumSides = " << fGeoEl->NSides() << " side = " << fSide << "\n";
371  }
372  if(!neighbour.Exists()) {
373  fGeoEl->SetSideDefined(fSide);
374  }
375 
376  TPZGeoElSide neighneigh, currentneigh;
377  neighneigh = neighbour.Neighbour();
378  currentneigh = Neighbour();
379  if (!neighneigh.Exists() && !currentneigh.Exists()) {
380  SetNeighbour(neighbour);
381  neighbour.SetNeighbour(*this);
382  } else if (neighneigh.Exists() && currentneigh.Exists()) {
383  // It would be convenient to check the consistency of both loops
384  // insert the connectivity between two independent loops
385  int a = NeighbourExists(neighbour);
386  int b = neighbour.NeighbourExists(*this);
387  if(
388  (a && !b) ||
389  (!a && b)
390  )
391  {
392  std::cout << "This element side : " << fSide << std::endl;
393  this->Element()->Print(std::cout);
394  cout << "\nNeighbour side :" << neighbour.Side() << std::endl;
395  neighbour.Element()->Print(std::cout);
396  PZError << "TPZGeoElSide::SetConnectivity Fourth untreated case, wrong data structure\n";
397  } else if(!a) {
398  SetNeighbour(neighneigh);
399  neighbour.SetNeighbour(currentneigh);
400  }
402  } else if (neighneigh.Exists()) {
403  // the neighbour already has a loop, insert this into the loop
404  SetNeighbour(neighneigh);
405  neighbour.SetNeighbour(*this);
407  } else if (currentneigh.Exists()) {
408  // this is already inserted in a loop, insert neighbour in the loop
409  SetNeighbour(neighbour);
410  neighbour.SetNeighbour(currentneigh);
411  }
412 }
413 
415 {
416  if(!fGeoEl) return;
417  TPZManVector<REAL,3> gelcenter(fGeoEl->Dimension());
418  fGeoEl->CenterPoint(fSide,gelcenter);
419  TPZTransform<> tr(Dimension(),fGeoEl->Dimension());
420  tr = fGeoEl->SideToSideTransform(fGeoEl->NSides()-1,fSide);
421  tr.Apply(gelcenter, center);
422 }
423 
426 {
427  if(!fGeoEl) return;
428  TPZManVector<REAL,3> gelcenter(fGeoEl->Dimension());
429  fGeoEl->CenterPoint(fSide,gelcenter);
430  fGeoEl->X(gelcenter, Xcenter);
431 }
432 
434  if(fSide < fGeoEl->NCornerNodes())
435  {
436  AllNeighbours(compneigh);
437  return;
438  }
439  int nsnodes = NSideNodes();
440  TPZStack<TPZGeoElSide> GeoElSideSet;
441  TPZStack<int> GeoElSet[27];
442  int in;
443  TPZManVector<int64_t> nodeindexes(nsnodes);
444  for(in=0; in<nsnodes; in++)
445  {
446  nodeindexes[in] = SideNodeIndex(in);
447  int locnod = fGeoEl->SideNodeLocIndex(fSide,in);
448  GeoElSideSet.Resize(0);
449  TPZGeoElSide locside(fGeoEl,locnod);
450  locside.AllNeighbours(GeoElSideSet);
451  int nel = GeoElSideSet.NElements();
452  int el;
453  for(el=0; el<nel; el++) {
454  GeoElSet[in].Push(GeoElSideSet[el].Element()->Index());
455  }
456  Sort<int>(GeoElSet[in]);
457  }
458  TPZStack<int> result;
459  switch(nsnodes) {
460  case 1:
461  {
462  result = GeoElSet[0];
463  }
464  break;
465  case 2:
466  Intersect<int,DEFAULTVEC_ALLOC>(GeoElSet[0],GeoElSet[1],result);
467  break;
468  case 3:
469  Intersect<int,DEFAULTVEC_ALLOC>(GeoElSet[0],GeoElSet[1],GeoElSet[2],result);
470  break;
471  case 4:
472  {
473  TPZStack<int> inter1, inter2;
474  Intersect<int,DEFAULTVEC_ALLOC>(GeoElSet[0],GeoElSet[2],inter1);
475  if(inter1.NElements()==0) break;
476  Intersect<int,DEFAULTVEC_ALLOC>(GeoElSet[1],GeoElSet[3],inter2);
477  if(inter2.NElements()==0) break;
478  Intersect<int,DEFAULTVEC_ALLOC>(inter1,inter2,result);
479  }
480  break;
481  default:
482  {
483  TPZStack<int> inter1, inter2;
484  inter1 = GeoElSet[0];
485  for(in=0; in<nsnodes-1; in++) {
486  inter2.Resize(0);
487  Intersect<int,DEFAULTVEC_ALLOC>(inter1,GeoElSet[in+1],inter2);
488  if(inter2.NElements() == 0) break;
489  inter1 = inter2;
490  }
491  result = inter2;
492  }
493  }
494  int el,nel = result.NElements();
495  TPZGeoMesh * geoMesh = fGeoEl->Mesh();
496  for(el=0; el<nel; el++) {
497  TPZGeoEl * gelResult = geoMesh->ElementVec()[result[el]];
498  int whichSd = gelResult->WhichSide(nodeindexes);
499  if(whichSd > 0)
500  {
501  compneigh.Push(TPZGeoElSide( gelResult, whichSd));
502  }
503  }
504 }
505 
506 
508 
509 #ifdef PZDEBUG
510  if(!NeighbourExists(neighbour))
511  {
512  stringstream sout;
513  sout << __PRETTY_FUNCTION__ << "Neighbour does not exist : expect trouble";
514  LOGPZ_ERROR(logger,sout.str());
515  TPZTransform<> toto;
516  return toto;
517  }
518 #endif
519  int sidedimension = Dimension();
520  TPZTransform<> tside(sidedimension);//transforma�o local
521  switch (sidedimension) {
522  case 0://canto para canto viz
523 
524  break;
525 
526  case 1://aresta para aresta viz
527  if (SideNodeIndex(0) == neighbour.SideNodeIndex(0))
528  {
529  tside.Mult()(0,0) = 1.;
530  }
531  else
532  {
533  tside.Mult()(0,0) = -1.;
534  }
535  break;
536 
537  case 2://transformacoes entre faces viz
538  int i;
539  //TPZCompEl *cel = Element()->Reference();
540  TPZVec<int> idto(0),idfrom(0);
541  if(Element()->NSideNodes(Side()) == 4) {//faces quadrilaterais
542  idto.Resize(4);
543  idfrom.Resize(4);
544  tside.Mult()(0,0) = 0;
545  tside.Mult()(1,1) = 0;
546  for(i=0;i<4;i++) idto[i]=neighbour.Element()->SideNodeIndex(neighbour.Side(),i);
547  for(i=0;i<4;i++) idfrom[i]=Element()->SideNodeIndex(Side(),i);
548 
549  int transid = Element()->GetTransformId2dQ(idfrom,idto);
550  tside.Mult()(0,0) = TPZShapeQuad::gTrans2dQ[transid][0][0];//cel->gTrans2dQ[transid][0][0];
551  tside.Mult()(0,1) = TPZShapeQuad::gTrans2dQ[transid][0][1];
552  tside.Mult()(1,0) = TPZShapeQuad::gTrans2dQ[transid][1][0];
553  tside.Mult()(1,1) = TPZShapeQuad::gTrans2dQ[transid][1][1];
554  } else if(Element()->NSideNodes(Side()) == 3) {//faces triangulares
555  idto.Resize(3);
556  idfrom.Resize(3);
557  tside.Mult()(0,0) = 0.;
558  tside.Mult()(1,1) = 0.;
559  for(i=0;i<3;i++) idfrom[i] = Element()->SideNodeIndex(Side(),i);
560  for(i=0;i<3;i++) idto[i] = neighbour.Element()->SideNodeIndex(neighbour.Side(),i);
561  int transid = Element()->GetTransformId2dT(idfrom,idto);
562  tside.Mult()(0,0) = TPZShapeTriang::gTrans2dT[transid][0][0];
563  tside.Mult()(0,1) = TPZShapeTriang::gTrans2dT[transid][0][1];
564  tside.Mult()(1,0) = TPZShapeTriang::gTrans2dT[transid][1][0];
565  tside.Mult()(1,1) = TPZShapeTriang::gTrans2dT[transid][1][1];
566  tside.Sum()(0,0) = TPZShapeTriang::gVet2dT[transid][0];
567  tside.Sum()(1,0) = TPZShapeTriang::gVet2dT[transid][1];
568  } else {
569  PZError << "TPZGeoElSide::NeighbourSideTransform : elemento desconhecido" << std::endl;
570  }
571  break;
572  }
573  return tside;
574 }
575 
577  if(gel == *this) return 1;
578  TPZGeoElSide neighbour = Neighbour();
579  if(!neighbour.Exists()) return 0;
580  while(neighbour != *this) {
581  if(gel == neighbour) return 1;
582  neighbour = neighbour.Neighbour();
583  }
584  return 0;
585 }
586 
588  if (!fGeoEl) return TPZCompElSide();
589  TPZCompElSide compside(fGeoEl->Reference(),fSide);
590  return compside;
591 }
592 
594  if (!fGeoEl) {
595  PZError << "TPZGeoElSide::Dimension : null element\n";
596  return -1;
597  }
598  return fGeoEl->SideDimension(fSide);
599 }
600 
602  //t : atual -> neighbour
603 #ifdef LOG4CXX
604  if (logger->isDebugEnabled()) {
605  std::stringstream sout;
606  sout << __FUNCTION__ << " this = \n";
607  Print(sout);
608  sout << "neighbour\n";
609  neighbour.Print(sout);
610  LOGPZ_DEBUG(logger, sout.str())
611  }
612 #endif
613  TPZGeoElSide father(*this);
614  if(!father.Exists()) {
615  PZError << "TPZGeoElSide::SideTransform3 I dont understand\n";
616  return;
617  }
618  while(father.Exists())
619  {
620  if(father.NeighbourExists(neighbour)) {
621  TPZTransform<> Temp = father.NeighbourSideTransform(neighbour);
622  // t = NeighbourSideTransform(neighbour).Multiply(t);
623  t = Temp.Multiply(t);
624  return;
625  } else {
626  TPZGeoElSide nextfather = father.StrictFather();
627  if(nextfather.Exists())
628  {
629  t = father.Element()->BuildTransform2(father.Side(),nextfather.Element(),t);
630  }
631  father = nextfather;
632  }
633  }
634  TPZGeoElSide start,neighbourwithfather,neighfather;
635  start = *this;
636  neighbourwithfather = *this;
637  do {
638  neighfather = neighbourwithfather.StrictFather();
639  if(!neighfather.Exists()) neighbourwithfather = neighbourwithfather.Neighbour();
640  } while(!neighfather.Exists() && neighbourwithfather.Exists() && neighbourwithfather != start);
641  int secondcase = 0;
642 
643  if(neighfather.Exists()) {
644  if(neighbourwithfather != start)
645  {
646  secondcase++;
647  t = start.NeighbourSideTransform(neighbourwithfather).Multiply(t);
648  }
649 #ifdef LOG4CXX
650  if (logger->isDebugEnabled()) {
651  std::stringstream sout;
652  sout << "neighbourwithfather\n";
653  neighbourwithfather.Print(sout);
654  sout << "neighbour\n";
655  neighbour.Print(sout);
656  LOGPZ_DEBUG(logger, sout.str())
657  }
658 #endif
659  neighbourwithfather.SideTransform3(neighbour,t);
660  return;
661  }
662  Print(std::cout);
663  neighbour.Print(std::cout);
664  fGeoEl->Print();
665  neighbour.Element()->Print();
666  PZError << "TPZGeoElSide:SideTranform3 did not find the neighbour\n";
667  return;
668 }
669 
671  int onlyinterpolated, int removeduplicates) {
672  if(!fGeoEl) return;
673  TPZCompElSide father = LowerLevelCompElementList2(onlyinterpolated);
674  if(father.Exists()) ellist.Push(father);
675  EqualLevelCompElementList(ellist,onlyinterpolated,removeduplicates);
676  HigherLevelCompElementList2(ellist,onlyinterpolated,removeduplicates);
677 }
678 
680  int onlyinterpolated, int removeduplicates) {
681 
682  TPZGeoElSide neighbour;
683  TPZCompElSide ref;
684  neighbour = Neighbour();
685  if(!neighbour.Exists()) return;
686 
687  while(neighbour.Element() != this->Element()) {
688  ref = neighbour.Reference();
689  if(ref.Element() && ref.Element() != Reference().Element() && (!onlyinterpolated || dynamic_cast<TPZInterpolatedElement*>(ref.Element()) )) {
690  elsidevec.Push(ref);
691  if(removeduplicates) return;
692  }
693  neighbour = neighbour.Neighbour();
694  }
695 }
696 
698  int onlymultiphysicelement, int removeduplicates) {
699 
700  TPZGeoElSide neighbour;
701  TPZCompElSide ref;
702  neighbour = Neighbour();
703  if(!neighbour.Exists()) return;
704 
705  while(neighbour.Element() != this->Element()) {
706  ref = neighbour.Reference();
707  if(ref.Element() && ref.Element() != Reference().Element() && (!onlymultiphysicelement || dynamic_cast<TPZMultiphysicsElement *>(ref.Element()) )) {
708  elsidevec.Push(ref);
709  if(removeduplicates) return;
710  }
711  neighbour = neighbour.Neighbour();
712  }
713  if(neighbour.Element() == this->Element())
714  {
715  ref = neighbour.Reference();
716  if(ref.Exists()) elsidevec.Push(ref);
717  }
718 }
719 
720 
722 
723  TPZStack<TPZGeoElSide> gelsides;
724  fGeoEl->AllHigherDimensionSides(fSide,2,gelsides);
725  int il,cap = gelsides.NElements();
726  for(il=0; il<cap; il++) {
727  TPZCompElSide cels = gelsides[il].Reference();
728  if(onlyinterpolated) {
729  TPZInterpolatedElement *cel = dynamic_cast<TPZInterpolatedElement *> (cels.Element());
730  if(!cel) continue;
731  //int64_t conind = cel->ConnectIndex(cels.Side());
732  int locconind = cel->MidSideConnectLocId(cels.Side());
733  int64_t conind = cel->ConnectIndex(locconind);
734  if(conind < 0) continue;
735  }
736  elsidevec.Push(cels);
737  }
738 
739 }
740 
741 int64_t TPZGeoElSide::Id() {
742  return fGeoEl->Id();
743 }
744 
746 void TPZGeoElSide::SetNeighbour(const TPZGeoElSide &neighbour) const {
747  if(!fGeoEl) {
748  PZError << "TPZGeoElSide::SetNeighbour. Don't exist geometrical element.\n";
749  return;
750  }
751  fGeoEl->SetNeighbour(fSide,neighbour);
752 }
753 
754 
756  if(fGeoEl != higherdimensionside.fGeoEl) {
757  PZError << "TPZGeoElSide::SideToSideTransform inconsistent id1 = " << fGeoEl->Id() <<
758  " id2 = " << higherdimensionside.fGeoEl->Id() << std::endl;
759  }
760  return fGeoEl->SideToSideTransform(fSide,higherdimensionside.fSide);
761 }
762 
764 {
765  TPZGeoEl * actGel = this->fGeoEl;
766  TPZGeoElSide side(*this);
767  while(actGel->Father())
768  {
769  side = actGel->Father2(side.Side());
770  actGel = actGel->Father();
771  }
772  return side;
773 }
774 
776 {
777 #ifdef PZDEBUG
778  PZError << __PRETTY_FUNCTION__ << "is deprecated. Use TPZGeoEl::YoungestChildren instead \n";
779  DebugStop();
780 #endif // PZDEBUG
781  if(this->Element()->HasSubElement() == false)
782  {
783  sonSides.Push(*this);
784  }
785  int dim = Dimension();
786  TPZStack<TPZGeoElSide> lowerSubelements;
787  fGeoEl->GetSubElements2(fSide,lowerSubelements,dim);
788  int nsub = lowerSubelements.size();
789  for (int s=0; s<nsub; s++) {
790  lowerSubelements[s].GetAllSiblings(sonSides);
791  }
792 }
793 
795 {
796  if(this->Element()->HasSubElement() == false)
797  {
798  sonSides.Push(*this);
799  }
800  int dim = Dimension();
801  TPZStack<TPZGeoElSide> lowerSubelements;
802  fGeoEl->GetSubElements2(fSide,lowerSubelements,dim);
803  int nsub = lowerSubelements.size();
804  for (int s=0; s<nsub; s++) {
805  lowerSubelements[s].YoungestChildren(sonSides);
806  }
807 }
808 
809 /*return 1 if the element has subelements along side*/
811  if(!fGeoEl) return 0;
812  return fGeoEl->HasSubElement();
813 }
814 /*return the number of nodes for a particular side*/
816  if(!fGeoEl) return 0;
817  return fGeoEl->NSideNodes(fSide);
818 }
819 
821 int64_t TPZGeoElSide::SideNodeIndex(int nodenum) const {
822  if(!fGeoEl) return -1;
823  return ( fGeoEl->SideNodeIndex(fSide,nodenum) );
824 }
825 
827 int64_t TPZGeoElSide::SideNodeLocIndex(int nodenum) const {
828  if(!fGeoEl) return -1;
829  return ( fGeoEl->SideNodeLocIndex(fSide,nodenum) );
830 }
831 
832 
834 {
835  // This method was modified to look for the father of any neighbouring element
836  // It is not sufficient to look for the father of the current element only, because a neighbour
837  // might have a father where the current element doesn t. This happens in the clone meshes. It probably
838  // will happen when working with interface elements or any situation where an element is inserted in an already refined mesh
839  TPZGeoElSide father,neighbour,start;
840  start = *this;
841  neighbour = Neighbour();
842  father = StrictFather();
843  int secondcase = 0;
844  while(!father.Exists() && neighbour.Exists() && neighbour != start) {
845  father = neighbour.StrictFather();
846  neighbour = neighbour.Neighbour();
847  secondcase++;
848  }
849  if(!father.Exists()) return TPZCompElSide();
851  //2003-12-03
852  if (father.Reference().Exists()) equal.Push(father.Reference());
853 
854  father.EqualLevelCompElementList(equal,onlyinterpolated,1);
855 
856 
857  while(father.Exists() && equal.NElements() == 0) {
858  neighbour = father.Neighbour();
859  start = father;
860  father = father.StrictFather();
861  while(!father.Exists() && neighbour.Exists() && neighbour != start) {
862  secondcase++;
863  father = neighbour.StrictFather();
864  neighbour = neighbour.Neighbour();
865  }
866  father.EqualLevelCompElementList(equal,onlyinterpolated,1);
867  //2003-12-03
868  if (father.Reference().Exists()) equal.Push(father.Reference());
869 
870  }
871  if(equal.NElements()) return equal[0];
872  return TPZCompElSide();
873 }
874 
876 {
877  if(!fGeoEl) return TPZGeoElSide();
878  return fGeoEl->Father2(fSide);
879 }
880 
882 {
883  TPZGeoElSide father = Father2();
884  int nfathsub = 0;
885  if(father.Exists()) nfathsub = father.fGeoEl->NSideSubElements(father.fSide);
886  while(father.Exists() && nfathsub == 1) {
887  father = father.Father2();
888  if(father.Exists()) nfathsub = father.fGeoEl->NSideSubElements(father.fSide);
889  }
890  return father;
891 }
892 
894 {
895  if(!fGeoEl || !fGeoEl->HasSubElement()) { // Jorge 10/01/2000
896 // comentei para poder acumular subelementos
897 // subelements.Resize(0);
898  return;
899  }
900  fGeoEl->GetSubElements2(fSide,subelements);
901 }
902 
903 void TPZGeoElSide::HigherLevelCompElementList2(TPZStack<TPZCompElSide> &elvec, int onlyinterpolated, int removeduplicates) {
904 
905  if(!Dimension()) return;
906  TPZGeoElSide neighbour(*this);
908  do {
909  if(neighbour.HasSubElement() && neighbour.NSubElements() > 1) {
910  neighbour.GetSubElements2(subel);
911  int nsub = subel.NElements();
912  int is;
913  for(is=0; is<nsub; is++) {
914  subel[is].EqualorHigherCompElementList2(elvec,onlyinterpolated,removeduplicates);
915  }
916  neighbour = *this;
917  } else {
918  neighbour = neighbour.Neighbour();
919  }
920  if(!neighbour.Exists()) break;
921  } while(neighbour != *this);
922 }
923 
924 void TPZGeoElSide::HigherLevelCompElementList3(TPZStack<TPZCompElSide> &elvec, int onlymultiphysicelement, int removeduplicates) {
925 
926  if(!Dimension()) return;
927  TPZGeoElSide neighbour(*this);
929  do {
930  if(neighbour.HasSubElement() && neighbour.NSubElements() > 1) {
931  neighbour.GetSubElements2(subel);
932  int nsub = subel.NElements();
933  int is;
934  for(is=0; is<nsub; is++) {
935  subel[is].EqualorHigherCompElementList3(elvec,onlymultiphysicelement,removeduplicates);
936  }
937  neighbour = *this;
938  } else {
939  neighbour = neighbour.Neighbour();
940  }
941  if(!neighbour.Exists()) break;
942  } while(neighbour != *this);
943 }
944 
945 
946 void TPZGeoElSide::EqualorHigherCompElementList2(TPZStack<TPZCompElSide> &celside, int onlyinterpolated, int removeduplicates){
947 
948 
949  int ncelsides = celside.NElements();
950  if(Reference().Exists()) {
951  celside.Push(Reference());
952  if(removeduplicates) {
953  return;
954  }
955  }
956  this->EqualLevelCompElementList(celside,onlyinterpolated,removeduplicates);
957  if(ncelsides != celside.NElements()) return;
958  TPZStack<TPZGeoElSide> gelsides;
959  TPZGeoElSide neighbour(*this);
960  do {
961  if(neighbour.HasSubElement() && neighbour.NSubElements() > 1) {
962  neighbour.GetSubElements2(gelsides);
963  int nsub = gelsides.NElements();
964  int is;
965  for(is=0; is<nsub; is++) {
966  gelsides[is].EqualorHigherCompElementList2(celside,onlyinterpolated,removeduplicates);
967  }
968  neighbour = *this;
969  } else {
970  neighbour = neighbour.Neighbour();
971  }
972  if(!neighbour.Exists()) break;
973  } while(neighbour != *this);
974 
975 }
976 
977 void TPZGeoElSide::EqualorHigherCompElementList3(TPZStack<TPZCompElSide> &celside, int onlymultiphysicelement, int removeduplicates){
978 
979 
980  int ncelsides = celside.NElements();
981  if(Reference().Exists()) {
982  celside.Push(Reference());
983  if(removeduplicates) {
984  return;
985  }
986  }
987  this->EqualLevelCompElementList3(celside,onlymultiphysicelement,removeduplicates);
988  if(ncelsides != celside.NElements()) return;
989  TPZStack<TPZGeoElSide> gelsides;
990  TPZGeoElSide neighbour(*this);
991  do {
992  if(neighbour.HasSubElement() && neighbour.NSubElements() > 1) {
993  neighbour.GetSubElements2(gelsides);
994  int nsub = gelsides.NElements();
995  int is;
996  for(is=0; is<nsub; is++) {
997  gelsides[is].EqualorHigherCompElementList3(celside,onlymultiphysicelement,removeduplicates);
998  }
999  neighbour = *this;
1000  } else {
1001  neighbour = neighbour.Neighbour();
1002  }
1003  if(!neighbour.Exists()) break;
1004  } while(neighbour != *this);
1005 
1006 }
1007 
1008 
1009 
1011 {
1012 
1013  if(!Exists()) return -1;
1014  return fGeoEl->NSideSubElements(fSide);
1015 }
1016 
1017 void TPZGeoElSide::BuildConnectivities(TPZVec<TPZGeoElSide> &sidevec,TPZVec<TPZGeoElSide> &neighvec){ //cout << "Sao iguais: acertar as vizinhancas!!!\n";
1022  int64_t size = sidevec.NElements();
1023  int64_t neighsize = neighvec.NElements();
1024  if(size!=neighsize || !size){
1025  PZError << "TPZGeoElSide::BuildConnectivities wrong vectors: abort!!!\n";
1026  DebugStop();
1027  }
1028  int64_t iv,ivn,side,neighside,sidedim,neighsidedim;
1029  TPZGeoElSide subside,neighsubside;
1030  for(iv=0;iv<size;iv++){
1031  subside = sidevec[iv];
1032  side = subside.Side();
1033  sidedim = subside.Dimension();
1034  TPZGeoEl *elside = subside.Element();
1035  for(ivn=0;ivn<neighsize;ivn++){
1036  neighsubside = neighvec[ivn];
1037  neighside = neighsubside.Side();
1038  neighsidedim = neighsubside.Dimension();
1039  TPZGeoEl *elneigh = neighsubside.Element();
1040  if(neighsidedim != sidedim) continue;
1041  //if(temp.Neighbour().Element()) continue;//?????? MELHORAR ESTA LINHA: pode ser viz. do irm�
1042  int in[4],face,i,j,num;
1043  int im[4],neighface;
1044  switch(sidedim){
1045  case 0://canto
1046  if(elside->SideNodeIndex(side,0) == elneigh->SideNodeIndex(neighside,0)){
1047  if(subside.NeighbourExists(neighsubside)) {
1048  cout << "TPZGeoElSide::BuildConnectivities the neighbour already exists?";
1049  } else {
1050  subside.SetConnectivity(neighsubside);
1051  }
1052  }
1053  break;
1054  case 1://aresta
1055  in[0] = elside->SideNodeIndex(side,0);
1056  in[1] = elside->SideNodeIndex(side,1);
1057  im[0] = elneigh->SideNodeIndex(neighside,0);
1058  im[1] = elneigh->SideNodeIndex(neighside,1);
1059  if( (in[0] == im[0] && in[1] == im[1]) || (in[0] == im[1] && in[1] == im[0]) ){
1060  if(subside.NeighbourExists(neighsubside)) {
1061  cout << "TPZGeoElSide::BuildConnectivities the neighbour already exists?";
1062  } else {
1063  subside.SetConnectivity(neighsubside);
1064  }
1065  }
1066  break;
1067  case 2://face
1068  //face = NSideNodes(side);//original substituida pela seguinte
1069  face = elside->NSideNodes(side);
1070  neighface = elneigh->NSideNodes(neighside);
1071  if(face!=neighface) break;
1072  for(i=0;i<face;i++){
1073  in[i] = elside->SideNodeIndex(side,i);
1074  im[i] = elneigh->SideNodeIndex(neighside,i);
1075  }
1076  if(face==3) in[3] = im[3] = -1;
1077  num = 0;
1078  for(i=0;i<4;i++) for(j=0;j<4;j++) if(in[i]==im[j]) num++;
1079  if(num==4){
1080  if(subside.NeighbourExists(neighsubside)) {
1081  cout << "TPZGeoElSide::BuildConnectivities the neighbour already exists?";
1082  } else {
1083  subside.SetConnectivity(neighsubside);
1084  }
1085  }
1086  break;
1087  default:
1088  PZError << "TPZGeoElSide::BuildConnectivities error !!!\n";
1089 
1090  }
1091  }
1092  }
1093 }
1094 
1095 std::ostream &operator << (std::ostream & out,const TPZGeoElSide &geoside){
1096  out << "TPZGeoElSide : side = " << geoside.Side() << std::endl ;
1097  if (geoside.Element()) geoside.Element()->Print(out);
1098  return out;
1099 }
1100 
1102 {
1103  if(!fGeoEl) return false;
1104  return fGeoEl->IsLinearMapping();
1105 }
1106 
1107 
1109 void TPZGeoElSide::Normal(TPZVec<REAL> &point, TPZGeoEl *LeftEl, TPZGeoEl *RightEl, TPZVec<REAL> &normal) const
1110 {
1111  normal.Resize(3,0.);
1112  normal.Fill(0.);
1113  int faceleft,faceright;
1114  int Leftdim = LeftEl->Dimension();
1115  int Rightdim = RightEl->Dimension();
1116 
1117  TPZManVector<REAL, 3> centleft(Leftdim),centright(Rightdim),result(3,0.),xint(3),xvolleft(3),xvolright(3),vec(3),rib(3);
1118  REAL normalize;
1119  int i;
1120 
1121  int InterfaceDimension = Dimension();
1122  TPZFNMatrix<9,REAL> axes(InterfaceDimension,3), jacobian(InterfaceDimension,InterfaceDimension),invjacobian(InterfaceDimension,InterfaceDimension);
1123  REAL detjac;
1124  this->Jacobian(point,jacobian,axes,detjac,invjacobian);
1125  faceleft = LeftEl->NSides()-1;//lado interior do elemento esquerdo
1126  faceright = RightEl->NSides()-1; // lado interior do element direito
1127  LeftEl->CenterPoint(faceleft,centleft);//ponto centro do elemento de volume
1128  RightEl->CenterPoint(faceright,centright);
1129  LeftEl->X(centleft,xvolleft);
1130  RightEl->X(centright,xvolright);
1131  for(i=0;i<3;i++) vec[i] = xvolright[i]-xvolleft[i];//nao deve ser nulo
1132 
1133 
1134  REAL vecnorm = sdot(vec, vec);
1135  if(vecnorm < 1.e-10)
1136  {
1137  LOGPZ_ERROR(logger,"Left and Right element centers coincide")
1138  vec[0]=1.;
1139  }
1140 
1141 
1142  switch(InterfaceDimension){
1143  case 0:
1144  normal[0] = vec[0];// a normal sempre aponta direcao positiva do eixo
1145  normal[1] = vec[1];
1146  normal[2] = vec[2];
1147 
1148  normalize = 0.;
1149  for(i=0;i<3;i++) normalize += normal[i]*normal[i];
1150  normalize = sqrt(normalize);
1151  if(!IsZero(normalize))
1152  for(i=0;i<3;i++) normal[i] = normal[i]/normalize;
1153 
1154  break;
1155  case 1:
1156  for(i=0;i<3;i++) rib[i] = axes(0,i);//direcao da aresta
1157  TPZNumeric::ProdVetorial(rib, vec, result);
1158  TPZNumeric::ProdVetorial(result, rib, normal);
1159  //normalizando a normal
1160  normalize = 0.;
1161  for(i=0;i<3;i++) normalize += normal[i]*normal[i];
1162  if(normalize == 0.0)
1163  {
1164  PZError << __PRETTY_FUNCTION__ << " null normal vetor\n";
1165 #ifdef LOG4CXX
1166  {
1167  std::stringstream sout;
1168  Print(sout);
1169  LOGPZ_DEBUG(logger,sout.str())
1170  }
1171 #endif
1172 
1173  DebugStop();
1174  }
1175  normalize = sqrt(normalize);
1176  for(i=0;i<3;i++) normal[i] = normal[i]/normalize;
1177  break;
1178  case 2:{
1179  TPZManVector<REAL,3> axes1(3), axes2(3);
1180  for(int iax = 0; iax < 3; iax++){
1181  axes1[iax] = axes(0,iax);
1182  axes2[iax] = axes(1,iax);
1183  }
1184  TPZNumeric::ProdVetorial(axes1,axes2,normal);
1185  }
1186  break;
1187  default:
1188  PZError << "TPZInterfaceElement::NormalToFace in case that not treated\n";
1189  normal.Resize(0);
1190  DebugStop();
1191  return;
1192  }
1193 
1194  //to guarantee the normal points from left to right neighbours:
1195  REAL dot = 0.;
1196  for(i=0; i<3; i++) dot += normal[i]*vec[i];
1197  if(dot < 0.) {
1198  for(i=0; i<3; i++) normal[i] = -normal[i];
1199  }
1200 
1201 }
1202 
1204 void TPZGeoElSide::Normal(TPZVec<REAL> &qsi_side, TPZVec<REAL> &normal) const{
1205 
1206  if (Dimension() != fGeoEl->Dimension()-1) {
1207  DebugStop();
1208  }
1209 
1210  normal.Resize(3,0.);
1211  normal.Fill(0.);
1212 
1213  TPZManVector<REAL, 3> center_dir(3,0.0), vol_center_x(3,0.0), x(3,0.0);
1214 
1215  TPZGeoElSide vol_side(fGeoEl,fGeoEl->NSides()-1);
1216  vol_side.CenterX(vol_center_x);
1217  this->X(qsi_side, x);
1218  center_dir = x - vol_center_x;
1219 
1220  int vol_dim = vol_side.Dimension();
1221  int side_dim = Dimension();
1222 
1223  TPZManVector<REAL, 3> qsi_vol(vol_dim,0.0);
1224  this->QsiElement(qsi_side, qsi_vol);
1225 
1226  TPZFMatrix<REAL> side_jac(side_dim,side_dim), side_inv_jac(side_dim,side_dim), side_axes(side_dim,3);
1227  TPZFMatrix<REAL> vol_jac(vol_dim,vol_dim), vol_inv_jac(vol_dim,vol_dim), vol_axes(vol_dim,3);
1228  REAL detjac;
1229 
1230  this->Jacobian(qsi_side, side_jac, side_axes, detjac, side_inv_jac);
1231  vol_side.Jacobian(qsi_vol, vol_jac, vol_axes, detjac, vol_inv_jac);
1232 
1233  switch (side_dim) {
1234  case 0:
1235  {
1236  for (unsigned int i = 0; i < 3; i++) {
1237  normal[i] = vol_axes(i,0);
1238  }
1239  }
1240  break;
1241  case 1:
1242  {
1243  TPZManVector<REAL, 3> v1(3,0.0),v2(3,0.0),v3(3,0.0),v4(3,0.0);
1244  for (unsigned int i = 0; i < 3 ; i++) {
1245  v1[i] = vol_axes(0,i);
1246  v2[i] = vol_axes(1,i);
1247  v3[i] = side_axes(0,i);
1248  }
1249 
1250  TPZNumeric::ProdVetorial(v1, v2, v4);
1251  TPZNumeric::ProdVetorial(v4, v3, normal);
1252 
1253  }
1254  break;
1255  case 2:
1256  {
1257  TPZManVector<REAL, 3> v1(3,0.0),v2(3,0.0);
1258  for (unsigned int i = 0; i < 3 ; i++) {
1259  v1[i] = side_axes(0,i);
1260  v2[i] = side_axes(1,i);
1261  }
1262  TPZNumeric::ProdVetorial(v1, v2, normal);
1263  }
1264  break;
1265  default:
1266  {
1267  DebugStop();
1268  }
1269  break;
1270  }
1271 
1273 
1274  //to guarantee the normal points from left to right neighbours:
1275  REAL dot = 0.;
1276  for(unsigned int i=0; i<3; i++) dot += normal[i]*center_dir[i];
1277  if(dot < 0.) {
1278  for(unsigned int i=0; i<3; i++) normal[i] = -normal[i];
1279  }
1280 
1281 }
1282 
1284 void TPZGeoElSide::Print(std::ostream &out) const
1285 {
1286  if(! fGeoEl)
1287  {
1288  out << "Null TPZGeoElSide\n";
1289  return;
1290  }
1291  out << "Element index " << fGeoEl->Index() << " Side " << fSide << " SideNode indexes " ;
1292  TPZManVector<REAL,3> center(Dimension(),0.), centerX(3,0.);
1293  CenterPoint(center);
1294  X(center, centerX);
1295  for (int i=0; i<NSideNodes(); i++) {
1296  out << SideNodeIndex(i) << " ";
1297  }
1298  out << "Center coordinate " << centerX << std::endl;
1299 }
1300 
1302  return this->Element()->CreateSideIntegrationRule(this->Side(), order);
1303 }
1304 
1305 int TPZGeoElSide::GelLocIndex(int index) const
1306 {
1307  if (!fGeoEl) {
1308  DebugStop();
1309  }
1310  return fGeoEl->SideNodeLocIndex(fSide,index);
1311 }
1312 
1314  return Hash("TPZGeoElSide");
1315 }
1316 
1317 void TPZGeoElSide::Read(TPZStream& buf, void* context) { //ok
1318  fGeoEl = dynamic_cast<TPZGeoEl*>(TPZPersistenceManager::GetInstance(&buf));
1319  buf.Read(&fSide);
1320 }
1321 
1322 void TPZGeoElSide::Write(TPZStream& buf, int withclassid) const { //ok
1324  buf.Write(&fSide);
1325 }
1326 
1328  if (this->fSide == -1 || this->fGeoElIndex == -1){
1329  return NULL;
1330  }
1331  return mesh->ElementVec()[this->fGeoElIndex];
1332 }
virtual int64_t SideNodeIndex(int side, int nodenum) const =0
Returns the index of the nodenum node of side.
TPZGeoEl * Element(const TPZGeoMesh *mesh) const
void EqualLevelCompElementList3(TPZStack< TPZCompElSide > &elsidevec, int onlymultiphysicelement, int removeduplicates)
Returns all connected computational elements which have equal level to the current element...
void EqualorHigherCompElementList2(TPZStack< TPZCompElSide > &celside, int onlyinterpolated, int removeduplicates)
Will return all elements of equal or higher level than than the current element.
void HigherLevelCompElementList2(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated, int removeduplicates)
Returns all connected computational elements which have level higher to the current element if onlyi...
int NSideNodes() const
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
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
virtual int NPoints() const =0
Returns number of points for the cubature rule related.
virtual void LowerDimensionSides(int side, TPZStack< int > &smallsides) const =0
Implements computational element and a side. Computational Element.
Definition: pzcompel.h:632
bool IsAncestor(TPZGeoEl *son, TPZGeoEl *father)
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
Contains declaration of TPZGeoElSide class which represents an element and its side, and TPZGeoElSideIndex class which represents an TPZGeoElSide index.
void GradX(TPZVec< REAL > &loc, TPZFMatrix< REAL > &gradx) const
X coordinate of a point loc of the side.
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.
int NSubElements()
return the number of element/side pairs which compose the current set of points
void SetElement(TPZGeoEl *geoel)
Definition: pzgeoelside.cpp:46
clarg::argInt nsub("-nsub", "number of substructs", 4)
int ClassId() const override
Define the class id associated with the class.
Definition: pzgeoelside.cpp:51
int GelLocIndex(int index) const
TPZGeoElSide StrictFather()
returns the father/side pair which contains this/side and is strictly larger than this/side ...
Contains declaration of TPZCompEl class which defines the interface of a computational element...
static REAL gTrans2dT[6][2][2]
Data structure which defines the triangle transformations.
TPZGeoElSide Father2() const
returns the father/side pair which contains the set of points associated with this/side ...
void CopyFrom(const TPZTransform< REAL > &cp)
Create a copy form a real transformation.
Definition: pztrnsform.cpp:62
groups all classes dedicated to the computation of shape functions
Definition: pzshapeextend.h:16
TPZGeoEl * fGeoEl
Definition: pzgeoelside.h:85
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
void GetSubElements2(TPZStack< TPZGeoElSide > &subelements)
build the list of element/side pairs which compose the current set of points
TPZCompElSide LowerLevelCompElementList2(int onlyinterpolated)
return the element/side pair which contains this/side and has a computational element associated ...
double sdot(TPZVec< T1 > &x, TPZVec< T1 > &y)
Performs a sdot operation: dot <- Transpose[x] * y.
Definition: pzvec_extras.h:98
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
Definition: fad.h:54
int HasSubElement()
Return 1 if the element has subelements along side.
TPZGeoElSide Neighbour() const
Definition: pzgeoel.h:754
virtual int NSides() const =0
Returns the number of connectivities of the element.
static TPZSavable * GetInstance(const int64_t &objId)
void CenterPoint(TPZVec< REAL > &center) const
return the coordinates of the center in master element space (associated with the side) ...
int ClassId() const override
Define the class id associated with the class.
TPZTransform< REAL > SideToSideTransform(TPZGeoElSide &higherdimensionside)
Compute the transformation between the master element space of one side of an element to the master e...
void SideTransform3(TPZGeoElSide neighbour, TPZTransform<> &t)
Accumulates the transformations from the current element/side to the neighbour/side.
void EqualLevelCompElementList(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated, int removeduplicates)
Returns all connected computational elements which have equal level to the current element...
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
int NNeighboursButThisElem(TPZGeoEl *thisElem)
Returns the number of neighbours, excluding the given element (thisElem)
static REAL gVet2dT[6][2]
Data structure which defines the triangle transformations.
virtual void YoungestChildren(TPZStack< TPZGeoElSide > &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.
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
It has the declaration of the TPZMultiphysicsCompEl class.
static void BuildConnectivities(TPZVec< TPZGeoElSide > &elvec, TPZVec< TPZGeoElSide > &neighvec)
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
void EqualorHigherCompElementList3(TPZStack< TPZCompElSide > &celside, int onlymultiphysicelement, int removeduplicates)
Will return all elements of equal or higher level than than the current element if onlymultiphysicel...
static void ProdVetorial(TPZVec< Tvar > &u, TPZVec< Tvar > &v, TPZVec< Tvar > &result)
Computes the vectorial product u x v.
Definition: pznumeric.cpp:96
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...
virtual int NSideSubElements(int side) const =0
Returns the number of subelements as returned by GetSubElements2(side)
bool IsLinearMapping() const
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZGeoElSide LowestFatherSide()
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
void AllNeighbours(TPZStack< TPZGeoElSide > &allneigh)
Returns the set of neighbours which can directly be accessed by the datastructure.
Definition: pzgeoelside.h:333
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
void Read(TPZStream &buf, void *context) override
read objects from the stream
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
int64_t Index() const
Returns the index of the element within the element vector of the mesh.
Definition: pzgeoel.h:730
TPZIntPoints * CreateIntegrationRule(int order)
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
void HigherDimensionElementList(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated)
Pushes all connected computational elements which have higher dimension than the current element/side...
void CenterX(TPZVec< REAL > &Xcenter) const
return the coordinates of the center of the side in real space
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
void X(TPZVec< REAL > &loc, TPZVec< REAL > &result) const
X coordinate of a point loc of the side.
virtual void GetAllSiblings(TPZStack< TPZGeoElSide > &unrefinedSons)
[deprecated] use YoungestChildren
TPZTransform< REAL > NeighbourSideTransform(const TPZGeoElSide &neighbour)
void RemoveConnectivity()
Remove the element from the connectivity loop.
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
const TPZFMatrix< T > & Sum() const
Definition: pztrnsform.h:66
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
int NSides() const
Returns the number of sides in which the current side can be decomposed.
void SetConnectivity(const TPZGeoElSide &neighbour) const
int Exists() const
Verifies if the object is non null (initialized)
Definition: pzcompel.h:697
TPZGeoElSideIndex()
Simple constructor.
Definition: pzgeoelside.h:356
static REAL gTrans2dQ[8][2][2]
Data structure which defines the quadrilateral transformations.
Definition: pzshapequad.h:180
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
virtual int NSideNodes(int side) const =0
Returns the number of nodes for a particular side.
A simple stack.
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
void ConnectedCompElementList(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated, int removeduplicates)
Returns all connected computational elements to the current element if onlyinterpolated == 1 only e...
virtual int64_t ConnectIndex(int i) const =0
Returns the index of the ith connectivity of the element.
int64_t Id()
TPZTransform< T > Multiply(TPZTransform< T > &right)
Multiply the transformation object (to the right) with right (Multiplying matrices) ...
Definition: pztrnsform.cpp:112
void ComputeNeighbours(TPZStack< TPZGeoElSide > &compneigh)
Returns the set of neighbours as computed by the intersection of neighbours along vertices...
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
int64_t SideNodeLocIndex(int nodenum) const
Returns the index of the local nodenum node of side.
void Normal(TPZVec< REAL > &point, TPZGeoEl *left, TPZGeoEl *right, TPZVec< REAL > &normal) const
compute the normal to the point from left to right neighbour
virtual int Dimension() const =0
Returns the dimension of the element.
TPZCompEl * Element() const
Gives a pointer to the reference computational element.
Definition: pzcompel.h:675
bool IsRelative(TPZGeoElSide other)
Checks whether other is a relative (son or ancestor) of this.
int64_t Id() const
Returns the Id of the element.
Definition: pzgeoel.h:223
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...
TPZCompElSide Reference() const
Returns a pointer to the elementside referenced by the geometric elementside.
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
REAL Area()
Area associated with the side.
int Dimension() const
the dimension associated with the element/side
bool ResetBlendConnectivity(const int64_t &index)
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.
void Print(std::ostream &out) const
print geometric characteristics of the element/side
bool IsAncestor(TPZGeoElSide other)
Checks whether other is an ancestor of this.
Definition: pzgeoelside.cpp:95
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
static int sidedimension[27]
Vector of the dimension for each side.
Definition: tpzcube.cpp:100
int Side() const
Definition: pzgeoelside.h:169
Contains declaration of TPZInterpolatedElement class which implements computational element of the in...
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
void SetNeighbour(const TPZGeoElSide &neighbour) const
Fill in the data structure for the neighbouring information.
int Side() const
Returns the side index.
Definition: pzcompel.h:688
int Exists() const
Definition: pzgeoelside.h:175
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzgeoelside.cpp:60
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
Contains the TPZTransform<> class which implements an affine transformation between points in paramet...
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
int NeighbourExists(const TPZGeoElSide &neighbour) const
Returns 1 if neighbour is a neighbour of the element along side.
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
void HigherLevelCompElementList3(TPZStack< TPZCompElSide > &elsidevec, int onlymultiphysicelement, int removeduplicates)
Returns all connected computational elements which have level higher to the current element if onlym...
void QsiElement(TPZVec< REAL > &qsi_side, TPZVec< REAL > &qsi_element) const
Parametric coordinate of a point loc of the side and return parametric element point.
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: pzgeoelside.cpp:55
Extra utilities for TPZVec. Implementations of the saxpy, sscal, sdot, intercept, max and min functio...
static void NormalizeVetor(TPZVec< Tvar > &vetor)
Normalizes the vector.
Definition: pznumeric.cpp:33
virtual void Print(std::ostream &out=std::cout)
Print all relevant data of the element to cout.
Definition: pzgeoel.cpp:238
Contains declaration of the TPZNumeric class which implements several methods to calculation.
TPZGeoEl * Father() const
Computes the normal vectors needed for forming HDiv vector valued shape functions.
Definition: pzgeoel.cpp:2574
Contains TPZShapeQuad class which implements the shape functions of a quadrilateral element...
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
int64_t SideNodeIndex(int nodenum) const
Returns the index of the nodenum node of side.
void InsertConnectivity(TPZGeoElSide &neighbour)
This method inserts the element/side and all lowerdimension sides into the connectivity loop...
static void WritePointer(const TPZSavable *obj, TPZStream *stream)
std::ostream & operator<<(std::ostream &out, const TPZGeoElSide &geoside)
Overload operator << to print geometric element side data.
Contains TPZShapeTriang class which implements the shape functions of a triangular element...
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
Implements computational element based on an interpolation space. Computational Element.
Definition: pzintel.h:27
virtual void Point(int i, TPZVec< REAL > &pos, REAL &w) const =0
Returns i-th point at master element and related weight.
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
virtual int MidSideConnectLocId(int is) const
Returns the local id of the connect in the middle of the side.
Definition: pzintel.cpp:83
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138