NeoPZ
tpzgeoblend.cpp
Go to the documentation of this file.
1 
5 #include "tpzgeoblend.h"
6 
7 #include "pzgeoelside.h"
8 #include "tpzgeoelmapped.h"
9 #include "pzgeoelrefless.h"
10 #include "tpzgeoelrefpattern.h"
11 #include "pzextractval.h"
12 #include "pzlog.h"
13 #include <sstream>
14 
15 #ifdef LOG4CXX
16 static LoggerPtr logger(Logger::getLogger("pz.mesh.geoblend"));
17 #endif
18 
19 template<class TGeo>
20 bool pzgeom::TPZGeoBlend<TGeo>::IsLinearMapping(int side) const
21 {
22  TPZStack<int> LowAllSides;
23  TGeo::LowerDimensionSides(side,LowAllSides);
24  if(side < 0 || side > TGeo::NSides-1)
25  {
26  DebugStop();
27  return 0;
28  }
29  bool straight = true;
30 
31  if(side >= TGeo::NNodes && fNeighbours[side-TGeo::NNodes].ElementIndex() != -1)
32  {
33  straight = false;
34  return straight;
35  }
36  for(int lowside = 0; lowside < LowAllSides.NElements(); lowside++)
37  {
38  if(LowAllSides[lowside] >= TGeo::NNodes && fNeighbours[LowAllSides[lowside]-TGeo::NNodes].ElementIndex() != -1)
39  {
40  straight = false;
41  return straight;
42  }
43  }
44  return straight;
45 }
46 
47 template <class TGeo>
48 bool pzgeom::TPZGeoBlend<TGeo>::ResetBlendConnectivity(const int64_t &side, const int64_t &index){
49  if(fNeighbours[side-TGeo::NNodes].ElementIndex() == index){
50  fNeighbours[side-TGeo::NNodes].SetElementIndex(-1);
51  return true;
52  }
53  return false;
54 }
55 
56 template <class TGeo>
57 void pzgeom::TPZGeoBlend<TGeo>::SetNeighbourInfo(int side, TPZGeoElSide &neigh, TPZTransform<> &trans)
58 {
59  if(!(fNeighbours[side-TGeo::NNodes].ElementIndex() != -1))
60  {
61  fNeighbours[side-TGeo::NNodes] = neigh;
62  fTrans[side - TGeo::NNodes] = trans;
63  }
64  else
65  {
66 #ifdef LOG4CXX
67  if(logger->isWarnEnabled())
68  {
69 
70  std::stringstream mess;
71  mess << "Trying to SetNeighbourInfo for an already set element\n";
72  mess << "* this * = " << __PRETTY_FUNCTION__ << "\n";
73  this->Print(mess);
74  mess << "* neigh * = \n";
75  neigh.Element()->Print(mess);
76  LOGPZ_WARN(logger,mess.str());
77  }
78 #endif
79  }
80 }
81 
82 template <class TGeo>
83 template<class T>
84 void pzgeom::TPZGeoBlend<TGeo>::GradX(TPZFMatrix<REAL> &coord, TPZVec<T> &xiInterior, TPZFMatrix<T> &gradx) const {
85 
86  TPZGeoEl &gel = *fGeoEl;
87  TPZGeoMesh *gmesh = gel.Mesh();
88 #ifdef LOG4CXX
89  const auto VAL_WIDTH = 10;
90  std::ostringstream soutLogDebug;
91  if(logger->isDebugEnabled())
92  { soutLogDebug<<"======================_______REF_1"<<std::endl;
93  soutLogDebug << "TPZGeoBlend<" <<MElementType_Name(TGeo::Type())<<">::GradX"<<std::endl;
94  soutLogDebug << "element id " <<gel.Id()<<std::endl;
95  soutLogDebug << "xi: ";
96  for(int i = 0; i < xiInterior.size(); i++) soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<xiInterior[i]<<"\t";
97  soutLogDebug<<std::endl;
98  }
99 #endif
100  const REAL zero = 1e-14;
101  gradx.Redim(3,TGeo::Dimension);
102 
103  if (fNeighbours[TGeo::NSides - 1 -TGeo::NNodes].ElementIndex() != -1){
104  TPZManVector<T, 3> neighXi;
105  TPZFNMatrix<9, T> dNeighXiDXi;
106  if (!MapToNeighSide(TGeo::NSides-1, TGeo::Dimension, xiInterior, neighXi, dNeighXiDXi)) {
107 #ifdef LOG4CXX2
108  if(logger->isDebugEnabled()) {
109  std::stringstream sout;
110  sout << "MapToNeighSide is singular for par " << xi << " and side " << TGeo::NSides-1 << ". Aborting...";
111  LOGPZ_DEBUG(logger,sout.str())
112  }
113 #endif
114  DebugStop();
115  }
116  TPZFNMatrix<9, T> gradNeigh;
117  Neighbour(TGeo::NSides-1, gmesh).GradX(neighXi, gradNeigh);
118  gradNeigh.Multiply(dNeighXiDXi,gradx);//gradNonLinSide = gradNeigh.dNeighXiDXi
119  return;
120  }
121 
127  TPZFNMatrix<9, T> phi(TGeo::NNodes, 1), dPhiDxi(TGeo::Dimension, TGeo::NNodes);
128  TGeo::TShape(xiInterior, phi, dPhiDxi);//gets the barycentric coordinates
129 
130 
131  TPZFNMatrix<45,T> gradXLin(3, TGeo::Dimension,(T)0);
132  for (int iNode = 0; iNode < TGeo::NNodes; iNode++) {//calculates the linear mapping
133  for(int x = 0; x < 3; x++) {
134  for (int xi = 0; xi < TGeo::Dimension; xi++) {
135  gradXLin(x, xi) += coord(x, iNode) * dPhiDxi(xi, iNode);
136  }
137  }
138  }
139  for(int i = 0; i < gradx.Rows(); i++){
140  for(int j = 0; j < gradx.Cols(); j++){
141  gradx(i,j) = gradXLin(i,j);
142  }
143  }
144 
145 #ifdef LOG4CXX
146  soutLogDebug<<"FIRST TERM"<<std::endl;
147  soutLogDebug << "gradient of linear mapping:\n";
148  if (logger->isDebugEnabled()) {
149  for (int i = 0; i < gradXLin.Rows(); i++) {
150  for (int j = 0; j < gradXLin.Cols(); j++) {
151  soutLogDebug << std::setw(VAL_WIDTH) << std::right << gradXLin(i, j);
152  if (j != gradXLin.Cols() - 1) soutLogDebug << "\t";
153  }
154  soutLogDebug << std::endl;
155  }
156  }
157 
158 #endif
159 
163  TPZVec<TPZFMatrix<T> > gradNonLinSideVec(TGeo::NSides - TGeo::NNodes,
164  TPZFNMatrix<27,T>(3, TGeo::Dimension, 0));
165  TPZVec<TPZFMatrix<T> > gradLinSideVec(TGeo::NSides - TGeo::NNodes,
166  TPZFNMatrix<27,T>(3, TGeo::Dimension, 0));
167  TPZManVector<T, 20> blendFactor(TGeo::NSides - TGeo::NNodes, (T)0);
168  TPZFNMatrix<27,T> dCorrFactorDxi(TGeo::NSides - TGeo::NNodes, TGeo::Dimension, (T) 0);
169  TPZFNMatrix<27, T> linearSideMappings(TGeo::NSides - TGeo::NNodes, 3, 0.);
170  TPZFNMatrix<27, T> nonLinearSideMappings(TGeo::NSides - TGeo::NNodes, 3, 0.);
171  for (int sideIndex = 0; sideIndex < TGeo::NSides - TGeo::NNodes - 1; sideIndex++) {
172  TPZManVector<T, 3> xiProjectedOverSide(TGeo::Dimension, 0);
173  TPZFMatrix<T> &gradNonLinSide = gradNonLinSideVec[sideIndex];
174  TPZFMatrix<T> &gradLinSide = gradLinSideVec[sideIndex];
175  gradLinSide.Redim(3,TGeo::Dimension);
176  gradNonLinSide.Redim(3,TGeo::Dimension);
177 
178  int side = TGeo::NNodes + sideIndex;
179  TPZGeoElSide gelside(fNeighbours[sideIndex], gmesh);
180  #ifdef LOG4CXX
181  if (logger->isDebugEnabled()) {
182  soutLogDebug << "================================" << std::endl;
183  soutLogDebug << "side: " << side << " is linear: ";
184  }
185  #endif
186  if (IsLinearMapping(side) || !gelside.Exists()) {
187  blendFactor[sideIndex] = 0;
188  #ifdef LOG4CXX
189  if (logger->isDebugEnabled()) {
190  if (IsLinearMapping(side)) soutLogDebug << "true" << std::endl;
191  else soutLogDebug << "false (no gelside) " << std::endl;
192  }
193  #endif
194  continue;
195  }
196  #ifdef LOG4CXX
197  if (logger->isDebugEnabled()) {
198  soutLogDebug << "false (gelside exists) " << std::endl;
199  }
200  #endif
201 
204  TPZManVector<T, 3> sideXi;
205  TPZFNMatrix<9, T> transfXiToSideXi;
206  bool regularMap = TGeo::CheckProjectionForSingularity(side, xiInterior);
207  if (!regularMap) {
208  #ifdef LOG4CXX
209  if (logger->isDebugEnabled()) {
210  soutLogDebug << "mapping is not regular. skipping side... ";
211 
212  }
213  #endif
214  continue;
215  }
216 
217  this->MapToSide(side, xiInterior, sideXi, transfXiToSideXi);
218  MElementType sideType = TGeo::Type(side);
219  const int nSideNodes = MElementType_NNodes(sideType);
220  const int sideDim = TGeo::SideDimension(side);
221  TPZFNMatrix<9, T> sidePhiVec(nSideNodes, 1),
222  dSidePhiDSideXiVec(sideDim, nSideNodes);
223  GetSideShapeFunction<TGeo>(side, sideXi, sidePhiVec, dSidePhiDSideXiVec);
224 
225  TPZFNMatrix<9, T> dXprojDxiSide(TGeo::Dimension,sideDim,(T)0),
226  dXiProjectedOverSideDxi(TGeo::Dimension,TGeo::Dimension,(T)0);
227  TPZManVector<REAL, 3> nodeCoord;
228  //calculation of transformation for the derivatives of sidephi
229 
230  TPZFNMatrix<9, T> gradLinSideXiSide(3,sideDim,(T)0);
231  for (int iNode = 0; iNode < nSideNodes; iNode++) {
232  const int currentNode = TGeo::SideNodeLocId(side, iNode);
233 
234  for (int x = 0; x < 3; x++) {
235  linearSideMappings(sideIndex, x) +=
236  coord(x, currentNode) * sidePhiVec(iNode, 0);
237  }
238 
239  TGeo::ParametricDomainNodeCoord(currentNode, nodeCoord);
240  for (int x = 0; x < TGeo::Dimension; x++) {
241  xiProjectedOverSide[x] += nodeCoord[x] * sidePhiVec(iNode, 0);
242  }
243 
244  for(int x = 0; x < 3; x++){
245  for(int xi = 0; xi < sideDim; xi++) {
246  gradLinSideXiSide(x,xi) += coord(x, currentNode) * dSidePhiDSideXiVec(xi,iNode);
247  }
248  }
249 
250  for (int i = 0; i < TGeo::Dimension; i++) {
251  for (int j = 0; j < sideDim; j++) {
252  dXprojDxiSide(i,j) += nodeCoord[i] * dSidePhiDSideXiVec(j,iNode);
253  }
254  }
255 
256  }
257 
258  gradLinSideXiSide.Multiply(transfXiToSideXi,gradLinSide);//gradLinSideTemp = gradLinSideXiSide.transfXiToSideXi
259  dXprojDxiSide.Multiply(transfXiToSideXi,dXiProjectedOverSideDxi);//dXprojDxiTemp = dXprojDxiSide.transfXiToSideXi
260 
261  #ifdef LOG4CXX
262  if (logger->isDebugEnabled()) {
263  soutLogDebug << "xi projection over side: ";
264  for (int x = 0; x < TGeo::Dimension; x++) soutLogDebug << xiProjectedOverSide[x] << "\t";
265  soutLogDebug << "\nGrad of projected point to side:\n";
266  for(int i = 0; i < dXiProjectedOverSideDxi.Rows(); i++){
267  for(int j = 0; j < dXiProjectedOverSideDxi.Cols(); j++){
268  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<dXiProjectedOverSideDxi(i,j);
269  if(j != dXiProjectedOverSideDxi.Cols() - 1) soutLogDebug<<"\t";
270  }
271  soutLogDebug<<std::endl;
272  }
273  soutLogDebug << std::endl << "xi projection in side domain:";
274  for (int x = 0; x < sideXi.size(); x++) soutLogDebug << sideXi[x] << "\t";
275  soutLogDebug << "\nTransformation from x side to xi interior:\n";
276  for(int i = 0; i < transfXiToSideXi.Rows(); i++){
277  for(int j = 0; j < transfXiToSideXi.Cols(); j++){
278  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<transfXiToSideXi(i,j);
279  if(j != transfXiToSideXi.Cols() - 1) soutLogDebug<<"\t";
280  }
281  soutLogDebug<<std::endl;
282  }
283  soutLogDebug << std::endl << "\nLinear mapping of projected point: ";
284  for (int x = 0; x < 3; x++) soutLogDebug << linearSideMappings(sideIndex, x) << "\t";
285  soutLogDebug << "\nGrad of linear mapping of point projected to side:\n";
286  for(int i = 0; i < gradLinSide.Rows(); i++){
287  for(int j = 0; j < gradLinSide.Cols(); j++){
288  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<gradLinSide(i,j);
289  if(j != gradLinSide.Cols() - 1) soutLogDebug<<"\t";
290  }
291  soutLogDebug<<std::endl;
292  }
293  }
294  #endif
295 
298  {
299  TPZManVector<T,3> dCorrFactor(TGeo::Dimension,(T)0);
300  TGeo::BlendFactorForSide(side, xiInterior, blendFactor[sideIndex], dCorrFactor);
301  for(int iXi = 0; iXi < TGeo::Dimension; iXi++){
302  dCorrFactorDxi(sideIndex,iXi) = dCorrFactor[iXi];
303  }
304  }
305 
306  int sidedim = gelside.Dimension();
307  TPZManVector<T, 3> neighXi;
308  TPZFNMatrix<9, T> dNeighXiDXi;
309  if (!MapToNeighSide(side, sidedim, xiInterior, neighXi, dNeighXiDXi)) {
310 #ifdef LOG4CXX2
311  if(logger->isDebugEnabled()) {
312  std::stringstream sout;
313  sout << "MapToNeighSide is singular for par " << par << " and side " << byside << " skipping the side ";
314  LOGPZ_DEBUG(logger,sout.str())
315  }
316 #endif
317  continue;
318  }
319  TPZManVector<T, 3> Xside(3, 0.);
320  Neighbour(side, gmesh).X(neighXi, Xside);
321  TPZFNMatrix<9, T> gradNeigh;
322  Neighbour(side, gmesh).GradX(neighXi,gradNeigh);
323  for (int x = 0; x < 3; x++) {
324  nonLinearSideMappings(sideIndex, x) = Xside[x];
325  }
326  gradNeigh.Multiply(dNeighXiDXi,gradNonLinSide);//gradNonLinSide = gradNeigh.dNeighXiDXi
327 
328 // else{
329  TPZStack<int> containedNodesInSide;
330  TGeo::LowerDimensionSides(side, containedNodesInSide, 0);
331 
332  TPZStack<int> allContainedSides;
333  TGeo::LowerDimensionSides(side, allContainedSides);
334  for (int subSideIndex = containedNodesInSide.NElements();
335  subSideIndex < allContainedSides.NElements(); subSideIndex++) {
336 // TPZFNMatrix<9, T> gradNonLinSubSide(3,TGeo::Dimension,(T)0);
337  const int subSide = allContainedSides[subSideIndex];
338 #ifdef LOG4CXX
339  if (logger->isDebugEnabled()) {
340  soutLogDebug << "\tSubside " << subSideIndex << " (global: " << subSide << ") is linear: ";
341  if (IsLinearMapping(subSide)) soutLogDebug << "true" << std::endl;
342  else soutLogDebug << "false" << std::endl;
343  }
344 #endif
345  if (IsLinearMapping(subSide)) continue;
346 
347  T blendFactorSide = -1;
348  TPZManVector<T,3> dCorrFactorSideDxiProj(TGeo::Dimension,(T)0);
349  TPZFNMatrix<3, T> dCorrFactorSideDxi(TGeo::Dimension,1,(T)0);
350  TGeo::BlendFactorForSide(subSide, xiProjectedOverSide, blendFactorSide, dCorrFactorSideDxiProj);
351 
352 #ifdef LOG4CXX
353  if (logger->isDebugEnabled()) {
354  soutLogDebug << "\tSubside influence :" << blendFactorSide << std::endl;
355  }
356 #endif
357  TPZFNMatrix<3, T> dCorrFactorSideDxiProjMat(1,TGeo::Dimension,(T)0);
358  for(int xi = 0; xi < TGeo::Dimension; xi++){
359  dCorrFactorSideDxiProjMat(0,xi) = dCorrFactorSideDxiProj[xi];
360  }
361 
362 // dCorrFactorSideDxi = dXiProjectedOverSideDxi. dCorrFactorSideDxiProjMat
363 // dXiProjectedOverSideDxi.Multiply(dCorrFactorSideDxiProjMat,dCorrFactorSideDxi);//mostrarphil
364 
365 // dCorrFactorSideDxi = dCorrFactorSideDxiProjMat . dXiProjectedOverSideDxi;
366  dCorrFactorSideDxiProjMat.Multiply(dXiProjectedOverSideDxi,dCorrFactorSideDxi);
367  #ifdef LOG4CXX
368  if (logger->isDebugEnabled()) {
369  soutLogDebug << "\n\tGrad of correction factor of point projected to side:\n";
370  for(int i = 0; i < dCorrFactorSideDxi.Rows(); i++){
371  soutLogDebug<<"\t";
372  for(int j = 0; j < dCorrFactorSideDxi.Cols(); j++){
373  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<dCorrFactorSideDxi(i,j);
374  if(j != dCorrFactorSideDxi.Cols() - 1) soutLogDebug<<"\t";
375  }
376  soutLogDebug<<std::endl;
377  }
378  }
379  #endif
380  TPZFMatrix<T> &gradNonLinSubSide = gradNonLinSideVec[subSide - TGeo::NNodes];
381  TPZFMatrix<T> &gradLinSubSide = gradLinSideVec[subSide - TGeo::NNodes];
382  for (int x = 0; x < 3; x++) {
383  nonLinearSideMappings(sideIndex, x) -=
384  blendFactorSide *
385  (nonLinearSideMappings(subSide - TGeo::NNodes, x) -
386  linearSideMappings(subSide - TGeo::NNodes, x));
387  for (int j = 0; j < TGeo::Dimension; j++) {
388  gradNonLinSide(x,j) -= blendFactorSide * (gradNonLinSubSide(x,j)-gradLinSubSide(x,j));
389  gradNonLinSide(x,j) -= dCorrFactorSideDxi(0,j) *
390  (nonLinearSideMappings(subSide - TGeo::NNodes, x) -
391  linearSideMappings(subSide - TGeo::NNodes, x));
392  }
393  }
394  }
395 // }
396 #ifdef LOG4CXX
397  if (logger->isDebugEnabled()) {
398  soutLogDebug << "Non-linear mapping: ";
399  for (int x = 0; x < 3; x++) soutLogDebug << nonLinearSideMappings(sideIndex, x) << "\t";
400  soutLogDebug << "\nGrad of non-linear mapping:\n";
401  for(int i = 0; i < gradNonLinSide.Rows(); i++){
402  for(int j = 0; j < gradNonLinSide.Cols(); j++){
403  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<gradNonLinSide(i,j);
404  if(j != gradNonLinSide.Cols() - 1) soutLogDebug<<"\t";
405  }
406  soutLogDebug<<std::endl;
407  }
408  soutLogDebug << std::endl;
409  soutLogDebug << "adding to result mapping of side: " << side << std::endl;
410  soutLogDebug << "\t\tcorrection factor: " << blendFactor[sideIndex] << std::endl;
411  soutLogDebug << "\t\tdCorrFactorDxi: ";
412  for(int i = 0; i < TGeo::Dimension; i++) soutLogDebug<<dCorrFactorDxi(sideIndex,i)<<"\t";
413  soutLogDebug << "\n";
414 
415  soutLogDebug<<"SECOND TERM (SIDE "<<sideIndex<<") :"<<std::endl;
416  soutLogDebug << "\t\tblendFactor * (gradNonLinSide-gradLinSide): " << std::endl;
417  for (int i = 0; i < 3; i++) {
418  for (int j = 0; j < TGeo::Dimension; j++) {
419  const T val = blendFactor[sideIndex] * (gradNonLinSide(i,j)-gradLinSide(i,j));
420  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<val;
421  if(j != TGeo::Dimension - 1) soutLogDebug<<"\t";
422  }
423  soutLogDebug << "\n";
424  }
425  soutLogDebug<<"THIRD TERM (SIDE "<<sideIndex<<") :"<<std::endl;
426  soutLogDebug << "\t\tdCorrFactorDxi *\n"
427  " (nonLinearSideMappings -\n"
428  " linearSideMappings): " << std::endl;
429  for (int i = 0; i < 3; i++) {
430  for (int j = 0; j < TGeo::Dimension; j++) {
431  const T val = dCorrFactorDxi(sideIndex,j) *
432  (nonLinearSideMappings(sideIndex, i) -
433  linearSideMappings(sideIndex,i));
434  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<val;
435  if(j != TGeo::Dimension - 1) soutLogDebug<<"\t";
436  }
437  soutLogDebug << "\n";
438  }
439 
440  soutLogDebug<<std::endl;
441  }
442 #endif
443 
444  for (int i = 0; i < 3; i++) {
445  for (int j = 0; j < TGeo::Dimension; j++) {
446  gradx(i,j) += blendFactor[sideIndex] * (gradNonLinSide(i,j)-gradLinSide(i,j));
447  gradx(i,j) += dCorrFactorDxi(sideIndex,j) *
448  (nonLinearSideMappings(sideIndex, i) -
449  linearSideMappings(sideIndex,i));
450  }
451  }
452  }
453 #ifdef LOG4CXX
454  if(logger->isDebugEnabled()){
455  soutLogDebug << "================================" <<std::endl;
456  soutLogDebug << "=============result=============" <<std::endl;
457  soutLogDebug <<"================================"<<std::endl;
458  for(int i = 0; i < gradx.Rows(); i++){
459  for(int j = 0; j < gradx.Cols(); j++){
460  soutLogDebug<<std::setw(VAL_WIDTH) << std::right<<gradx(i,j);
461  if(j != gradx.Cols() - 1) soutLogDebug<<"\t";
462  }
463  soutLogDebug<<std::endl;
464  }
465  soutLogDebug << "================================" <<std::endl;
466  soutLogDebug << "================================" <<std::endl;
467  soutLogDebug << "================================" <<std::endl;
468  LOGPZ_DEBUG(logger,soutLogDebug.str())
469  soutLogDebug.str("");
470  }
471 #endif
472 #ifdef PZDEBUG
473  for(int i = 0; i < gradx.Rows();i++){
474  for(int j = 0; j < gradx.Cols(); j++) {
475  if (std::isnan(TPZExtractVal::val(gradx(i,j)))) {
476  DebugStop();
477  }
478  }
479  }
480 #endif
481 
482 }
483 
484 template<class TGeo>
485 template<class T>
486 void pzgeom::TPZGeoBlend<TGeo>::X(TPZFMatrix<REAL> &coord, TPZVec<T> &xi, TPZVec<T> &result) const {
487  TPZGeoEl &gel = *fGeoEl;
488  TPZGeoMesh *gmesh = gel.Mesh();
489  TPZManVector<T,3> notUsedHereVec(3,(T)0);
490  TPZFNMatrix<9,T> notUsedHereMat(TGeo::NSides, TGeo::NSides,(T)0);//since some methods dont resize the matrix, it is
491  // bigger than needed.
492 
493  #ifdef LOG4CXX
494  std::ostringstream soutLogDebug;
495  if(logger->isDebugEnabled())
496  {
497  soutLogDebug << "TPZGeoBlend<" <<MElementType_Name(TGeo::Type())<<">::X2_______REF_1"<<std::endl;
498  soutLogDebug << "element id " <<gel.Id()<<std::endl;
499  soutLogDebug << "xi: ";
500  for(int i = 0; i < xi.size(); i++) soutLogDebug<<xi[i]<<"\n";
501  soutLogDebug<<std::endl;
502  }
503  #endif
504  const REAL zero = 1e-14;
505  result.Resize(3);
506  result.Fill(0);
507 
508  if (fNeighbours[TGeo::NSides - 1 -TGeo::NNodes].ElementIndex() != -1){
509  TPZManVector<T, 3> neighXi;
510  if (!MapToNeighSide(TGeo::NSides-1, TGeo::Dimension, xi, neighXi, notUsedHereMat)) {
511 #ifdef LOG4CXX2
512  if(logger->isDebugEnabled()) {
513  std::stringstream sout;
514  sout << "MapToNeighSide is singular for par " << xi << " and side " << TGeo::NSides-1 << ". Aborting...";
515  LOGPZ_DEBUG(logger,sout.str())
516  }
517 #endif
518  DebugStop();
519  }
520  Neighbour(TGeo::NSides-1, gmesh).X(neighXi, result);
521  return;
522  }
523 
529  TPZFNMatrix<9, T> phi(TGeo::NNodes, 1);
530  TGeo::TShape(xi, phi, notUsedHereMat);//gets the barycentric coordinates
531 
532 
533  for (int iNode = 0; iNode < TGeo::NNodes; iNode++) {//calculates the linear mapping
534  for (int x = 0; x < 3; x++) {
535  result[x] += coord(x, iNode) * phi(iNode, 0);
536  }
537  }
538  #ifdef LOG4CXX
539  if(logger->isDebugEnabled())
540  {
541  soutLogDebug << "Linear mapping:\n"<<std::endl;
542  for(int i = 0; i < result.size(); i++) soutLogDebug<<result[i]<<"\n";
543  soutLogDebug<<std::endl;
544  }
545  #endif
546 
550  TPZManVector<bool, 20> isRegularMapping(TGeo::NSides - TGeo::NNodes, 0.);
551  for (int sideIndex = 0; sideIndex < TGeo::NSides - TGeo::NNodes - 1; sideIndex++) {
552  int side = TGeo::NNodes + sideIndex;
553  isRegularMapping[sideIndex] = TGeo::CheckProjectionForSingularity(side,xi);
554  }
555 
556  TPZManVector<T, 20> blendFactor(TGeo::NSides - TGeo::NNodes, 0.);
557  TPZFNMatrix<27, T> projectedPointOverSide(TGeo::NSides - TGeo::NNodes, TGeo::Dimension, 0.);
558  TPZFNMatrix<27, T> linearSideMappings(TGeo::NSides - TGeo::NNodes, 3, 0.);
559  TPZFNMatrix<27, T> nonLinearSideMappings(TGeo::NSides - TGeo::NNodes, 3, 0.);
560  for (int sideIndex = 0; sideIndex < TGeo::NSides - TGeo::NNodes - 1; sideIndex++) {
561  int side = TGeo::NNodes + sideIndex;
562  if(!isRegularMapping[sideIndex]) {
563  #ifdef LOG4CXX
564  if(logger->isDebugEnabled()){
565  soutLogDebug <<"mapping is not regular. skipping side... ";
566 
567  }
568  #endif
569  continue;
570  }
571  TPZGeoElSide gelside(fNeighbours[sideIndex], gmesh);
572  #ifdef LOG4CXX
573  if(logger->isDebugEnabled())
574  {
575  soutLogDebug << "================================" <<std::endl;
576  soutLogDebug << "side: "<<side<<" is linear: ";
577  }
578  #endif
579  if (IsLinearMapping(side) || !gelside.Exists()) {
580  blendFactor[sideIndex] = 0;
581  #ifdef LOG4CXX
582  if(logger->isDebugEnabled()){
583  if( IsLinearMapping(side) ) soutLogDebug <<"true"<<std::endl;
584  else soutLogDebug <<"false (no gelside) "<<std::endl;
585  }
586  #endif
587  continue;
588  }
589  #ifdef LOG4CXX
590  if(logger->isDebugEnabled()){
591  soutLogDebug <<"false (gelside exists) "<<std::endl;
592  }
593  #endif
594 
597  TPZManVector<T, 3> sideXi;
598  this->MapToSide(side, xi, sideXi, notUsedHereMat);
599 
600  MElementType sideType = TGeo::Type(side);
601  const int nSideNodes = MElementType_NNodes(sideType);
602  TPZFNMatrix<9, T> sidePhi(nSideNodes, 1);
603  notUsedHereMat.Resize(TGeo::Dimension,nSideNodes);
604  GetSideShapeFunction<TGeo>(side, sideXi, sidePhi, notUsedHereMat);
605 
606  TPZManVector<REAL,3> nodeCoord;
607  for (int iNode = 0; iNode < nSideNodes; iNode++) {
608  const int currentNode = TGeo::SideNodeLocId(side, iNode);
609 
610  for (int x = 0; x < 3; x++) {
611  linearSideMappings(sideIndex, x) +=
612  coord(x, currentNode) * sidePhi(iNode, 0);
613  }
614  TGeo::ParametricDomainNodeCoord(currentNode,nodeCoord);
615  for (int x = 0; x < TGeo::Dimension; x++) {
616  projectedPointOverSide(sideIndex,x) += nodeCoord[x] * sidePhi(iNode,0);
617  }
618  }
622  #ifdef LOG4CXX
623  if(logger->isDebugEnabled()){
624  soutLogDebug <<"xi projection over side:\n";
625  for (int x = 0; x < TGeo::Dimension; x++) soutLogDebug<<projectedPointOverSide(sideIndex,x)<<"\n";
626  soutLogDebug<<std::endl<<"xi projection in side domain:\n";
627  for (int x = 0; x < sideXi.size(); x++) soutLogDebug<<sideXi[x]<<"\n";
628  soutLogDebug<<std::endl<<"linear mapping of projected point:\n";
629  for (int x = 0; x < 3; x++) soutLogDebug<<linearSideMappings(sideIndex, x)<<"\n";
630  }
631  #endif
632 
633  TGeo::BlendFactorForSide(side, xi, blendFactor[sideIndex], notUsedHereVec);
634  int sidedim = gelside.Dimension();
635  TPZManVector<T, 3> neighXi;
636  if (!MapToNeighSide(side, sidedim, xi, neighXi, notUsedHereMat)) {
637 #ifdef LOG4CXX2
638  if(logger->isDebugEnabled()) {
639  std::stringstream sout;
640  sout << "MapToNeighSide is singular for par " << par << " and side " << byside << " skipping the side ";
641  LOGPZ_DEBUG(logger,sout.str())
642  }
643 #endif
644  continue;
645  }
646  TPZManVector<T, 3> Xside(3, 0.);
647  Neighbour(side, gmesh).X(neighXi, Xside);
648  for (int x = 0; x < 3; x++) {
649  nonLinearSideMappings(sideIndex, x) = Xside[x];
650  }
651 
652 // else{
653  TPZStack<int> containedNodesInSide;
654  TGeo::LowerDimensionSides(side, containedNodesInSide, 0);
655 
656  TPZStack<int> allContainedSides;
657  TGeo::LowerDimensionSides(side, allContainedSides);
658  for (int subSideIndex = containedNodesInSide.NElements();
659  subSideIndex < allContainedSides.NElements(); subSideIndex++) {
660  const int subSide = allContainedSides[subSideIndex];
661 #ifdef LOG4CXX
662  if (logger->isDebugEnabled()) {
663  soutLogDebug << "\tSubside " << subSideIndex << " (global: " << subSide << ") is linear: ";
664  if (IsLinearMapping(subSide)) soutLogDebug << "true" << std::endl;
665  else soutLogDebug << "false" << std::endl;
666  }
667 #endif
668  if(!isRegularMapping[subSide - TGeo::NNodes]) {
669  #ifdef LOG4CXX
670  if(logger->isDebugEnabled()){
671  soutLogDebug <<"mapping of subside is not regular. skipping side... ";
672 
673  }
674  #endif
675  continue;
676  }
677  if (IsLinearMapping(subSide)) continue;
678  TPZManVector<T, 3> projectedPoint(TGeo::Dimension, -1);
679  for (int x = 0; x < TGeo::Dimension; x++) {
680  projectedPoint[x] = projectedPointOverSide(sideIndex, x);
681  }
682 
683  T blendFactorSide = -1;
684  TGeo::BlendFactorForSide(subSide, projectedPoint, blendFactorSide,notUsedHereVec);
685 
686 #ifdef LOG4CXX
687  if (logger->isDebugEnabled()) {
688  soutLogDebug << "\tSubside influence :" << blendFactorSide << std::endl;
689  }
690 #endif
691  for (int x = 0; x < 3; x++) {
692  nonLinearSideMappings(sideIndex, x) -=
693  blendFactorSide *
694  (nonLinearSideMappings(subSide - TGeo::NNodes, x) -
695  linearSideMappings(subSide - TGeo::NNodes, x));
696  }
697 
698  }
699 // }
700 #ifdef LOG4CXX
701  if (logger->isDebugEnabled()) {
702  soutLogDebug << "Non-linear mapping\n";
703  for (int x = 0; x < 3; x++) soutLogDebug << nonLinearSideMappings(sideIndex, x) << "\n";
704  soutLogDebug << std::endl;
705 
706  soutLogDebug << "adding to result mapping of side: " << side << std::endl;
707  soutLogDebug << "\t\tblend factor: " << blendFactor[sideIndex] << std::endl;
708 // LOGPZ_DEBUG(logger,soutLogDebug.str())
709 // soutLogDebug.str("");
710  }
711 #endif
712  for (int x = 0; x < 3; x++) {
713  result[x] += blendFactor[sideIndex] *
714  (nonLinearSideMappings(sideIndex, x) - linearSideMappings(sideIndex, x));
715  }
716  }
717 #ifdef LOG4CXX
718  if(logger->isDebugEnabled()){
719  soutLogDebug << "================================" <<std::endl;
720  soutLogDebug << "=============result=============" <<std::endl;
721  soutLogDebug <<"================================"<<std::endl;
722  for (int x = 0; x < 3; x++) soutLogDebug<<result[x]<<"\n";
723  soutLogDebug << "\n================================" <<std::endl;
724  soutLogDebug << "================================" <<std::endl;
725  soutLogDebug << "================================" <<std::endl;
726  LOGPZ_DEBUG(logger,soutLogDebug.str())
727  soutLogDebug.str("");
728  }
729 #endif
730 }
731 
732 template <class TGeo>
733 void pzgeom::TPZGeoBlend<TGeo>::Jacobian(TPZFMatrix<REAL> &coord, TPZVec<REAL>& par, TPZFMatrix<REAL> &jacobian, TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv) const
734 {
735 
736  TPZGeoEl &gel = *fGeoEl;
737  TPZManVector<REAL,3> NeighPar, SidePar, Xside(3,0.), XNode(3,0.);
738  int majorSide = TGeo::NSides - 1;
739 
740  TPZManVector<REAL> SidesCounter(TGeo::NSides,0);
741  TPZStack<int> LowNodeSides, LowAllSides;
742 
743 #ifdef LOG4CXX
744  if(logger->isDebugEnabled())
745  {
746  std::stringstream sout;
747  sout << "input parameter par " << par;
748  LOGPZ_DEBUG(logger,sout.str())
749  }
750 #endif
751 
752  TPZFNMatrix<24> blend(TGeo::NNodes,1), Dblend(TGeo::Dimension,TGeo::NNodes);
753  TGeo::Shape(par,blend,Dblend);
754 
755  TPZFNMatrix<9> J1, J2, Ax, JacTemp(3,TGeo::Dimension, 0.), Jneighbourhood;
756  REAL Det;
757  TPZGeoMesh *gmesh = gel.Mesh();
758  for(int byside = majorSide; byside >= TGeo::NNodes; byside--)
759  {
760  TPZGeoElSide neighbyside = Neighbour(byside, gmesh);
761  if(neighbyside.Exists())
762  {
763  TGeo::LowerDimensionSides(byside,LowNodeSides,0);
764  TGeo::LowerDimensionSides(byside,LowAllSides);
765  int dim = Neighbour(byside, gmesh).Dimension();
766  TPZFNMatrix<9> Inv(dim,dim);
767  int sidedim = neighbyside.Dimension();
768  if(!MapToNeighSide(byside,sidedim,par,NeighPar, Jneighbourhood))
769  {
770  continue;
771  }
772  Neighbour(byside,gmesh).X(NeighPar,Xside);
773  Neighbour(byside,gmesh).Jacobian(NeighPar,J1,Ax,Det,Inv);
774  Ax.Transpose();
775  Ax.Multiply(J1,J2);
776 
777 #ifdef LOG4CXX
778  if(logger->isDebugEnabled())
779  {
780  std::stringstream sout;
781  sout << "byside " << byside << std::endl;
782  sout << "side of the neighbour " << Neighbour(byside,gmesh).Side() << std::endl;
783  Neighbour(byside,gmesh).Element()->Print(sout);
784  sout << "neighbour parameter(NeighPar) " << NeighPar << std::endl;
785  sout << "Jacobian of the map(Jneighborhood) " << Jneighbourhood << std::endl;
786  sout << "Xside " << Xside << std::endl;
787  sout << "jacobian neighbour(J1) " << J1 << std::endl;
788  Ax.Print("Ax of the neighbour (Ax) ",sout);
789  sout << "jacobian of the neighbour multiplied by the axes(J2) " << J2 << std::endl;
790  LOGPZ_DEBUG(logger,sout.str())
791  }
792 #endif
793 
794  J2.Multiply(Jneighbourhood,J1);
795 
796 #ifdef LOG4CXX
797  if(logger->isDebugEnabled())
798  {
799  std::stringstream sout;
800  sout << "acumulated jacobian(J1) " << J1 << std::endl;
801  LOGPZ_DEBUG(logger,sout.str())
802  }
803 #endif
804 
805  REAL blendTemp = 0.;
806  TPZManVector<REAL,3> DblendTemp(TGeo::Dimension,0.);
807  for(int a = 0; a < LowNodeSides.NElements(); a++)
808  {
809  TPZManVector<REAL> parChanged(par.NElements());
810  TGeo::Shape(parChanged,blend,Dblend);
811 
812  blendTemp += blend(LowNodeSides[a],0);
813  for(int b = 0; b < TGeo::Dimension; b++)
814  {
815  DblendTemp[b] += Dblend(b,LowNodeSides[a]);
816  }
817  }
818 
819  for(int a = 0; a < 3; a++)
820  {
821  for(int b = 0; b < TGeo::Dimension; b++)
822  {
823  JacTemp(a,b) += (1 - SidesCounter[byside]) * (J1(a,b)*blendTemp + Xside[a]*DblendTemp[b]);
824  }
825  }
826  for(int a = 0; a < LowAllSides.NElements(); a++)
827  {
828  SidesCounter[LowAllSides[a]] += (1 - SidesCounter[byside]);
829  }
830  }
831  }
832 
833 #ifdef LOG4CXX
834  if(logger->isDebugEnabled())
835  {
836  std::stringstream sout;
837  JacTemp.Print("Jabobian before contributing the nodes",sout);
838  sout << "SidesCounter " << SidesCounter << std::endl;
839  sout << "DBlend " << Dblend << std::endl;
840  sout << "NodeCoord " << coord << std::endl;
841  LOGPZ_DEBUG(logger,sout.str())
842  }
843 #endif
844 
845  for(int a = 0; a < TGeo::NNodes; a++)
846  {
847  for(int b = 0; b < 3; b++)
848  {
849  for(int c = 0; c < TGeo::Dimension; c++)
850  {
851  JacTemp(b,c) += (1 - SidesCounter[a]) * coord(b,a)*Dblend(c,a);
852  }
853  }
854  }
855 
856  if(TGeo::Dimension == 1)
857  {
858  TPZFNMatrix<9> axest;
859  JacTemp.GramSchmidt(axest,jacobian);
860  axest.Transpose(&axes);
861  jacinv.Resize(jacobian.Rows(), jacobian.Cols());
862  detjac = jacobian(0,0);
863  if(IsZero(detjac)){
864  detjac = ZeroTolerance();
865  }
866  jacinv(0,0) = 1./detjac;
867  } else if(TGeo::Dimension == 2)
868  {
869  TPZFNMatrix<9> axest;
870  JacTemp.GramSchmidt(axest,jacobian);
871  axest.Transpose(&axes);
872 
873  jacinv.Resize(jacobian.Rows(), jacobian.Cols());
874 
875  detjac = jacobian(0,0)*jacobian(1,1) - jacobian(1,0)*jacobian(0,1);
876  if(IsZero(detjac)){
877  detjac = ZeroTolerance();
878  }
879  jacinv(0,0) = jacobian(1,1) / detjac;
880  jacinv(1,1) = jacobian(0,0) / detjac;
881  jacinv(0,1) = -jacobian(0,1) / detjac;
882  jacinv(1,0) = -jacobian(1,0) / detjac;
883  }
884  else
885  {
886  jacobian = JacTemp;
887  axes.Resize(3,3); axes.Zero();
888 
889  jacinv.Resize(jacobian.Rows(), jacobian.Cols());
890 
891  axes(0,0) = 1.; axes(1,1) = 1.; axes(2,2) = 1.;
892  detjac = -jacobian(0,2)*jacobian(1,1)*jacobian(2,0);//- a02 a11 a20
893  detjac += jacobian(0,1)*jacobian(1,2)*jacobian(2,0);//+ a01 a12 a20
894  detjac += jacobian(0,2)*jacobian(1,0)*jacobian(2,1);//+ a02 a10 a21
895  detjac -= jacobian(0,0)*jacobian(1,2)*jacobian(2,1);//- a00 a12 a21
896  detjac -= jacobian(0,1)*jacobian(1,0)*jacobian(2,2);//- a01 a10 a22
897  detjac += jacobian(0,0)*jacobian(1,1)*jacobian(2,2);//+ a00 a11 a22
898 
899  if(IsZero(detjac))
900  {
901 #ifdef LOG4CXX
902  if(logger->isDebugEnabled())
903  {
904  std::stringstream sout;
905  sout << "Singular Jacobian " << detjac;
906  LOGPZ_ERROR(logger, sout.str())
907  }
908 #endif
909  detjac = ZeroTolerance();
910  }
911 
912  jacinv(0,0) = (-jacobian(1,2)*jacobian(2,1)+jacobian(1,1)*jacobian(2,2)) / detjac;//-a12 a21 + a11 a22
913  jacinv(0,1) = ( jacobian(0,2)*jacobian(2,1)-jacobian(0,1)*jacobian(2,2)) / detjac;// a02 a21 - a01 a22
914  jacinv(0,2) = (-jacobian(0,2)*jacobian(1,1)+jacobian(0,1)*jacobian(1,2)) / detjac;//-a02 a11 + a01 a12
915  jacinv(1,0) = ( jacobian(1,2)*jacobian(2,0)-jacobian(1,0)*jacobian(2,2)) / detjac;// a12 a20 - a10 a22
916  jacinv(1,1) = (-jacobian(0,2)*jacobian(2,0)+jacobian(0,0)*jacobian(2,2)) / detjac;//-a02 a20 + a00 a22
917  jacinv(1,2) = ( jacobian(0,2)*jacobian(1,0)-jacobian(0,0)*jacobian(1,2)) / detjac;// a02 a10 - a00 a12
918  jacinv(2,0) = (-jacobian(1,1)*jacobian(2,0)+jacobian(1,0)*jacobian(2,1)) / detjac;//-a11 a20 + a10 a21
919  jacinv(2,1) = ( jacobian(0,1)*jacobian(2,0)-jacobian(0,0)*jacobian(2,1)) / detjac;// a01 a20 - a00 a21
920  jacinv(2,2) = (-jacobian(0,1)*jacobian(1,0)+jacobian(0,0)*jacobian(1,1)) / detjac;//-a01 a10 + a00 a11
921  }
922 }
923 
924 template <class TGeo>
925 void pzgeom::TPZGeoBlend<TGeo>::Print(std::ostream &out) const
926 {
927  TGeo::Print(out);
928  out << "Neighbours/transformations used for mapping the sides :\n";
929  int is;
930  for(is=TGeo::NNodes; is<TGeo::NSides; is++)
931  {
932  out << "Side: " << is << " El/side: " << fNeighbours[is-TGeo::NNodes].ElementIndex() << ":" <<
933  fNeighbours[is-TGeo::NNodes].Side() << '\n';
934  }
935 }
936 
937 template <class TGeo>
938 void pzgeom::TPZGeoBlend<TGeo>::Initialize(TPZGeoEl *refel)
939 {
940  fGeoEl = refel;
941  for(int byside = TGeo::NNodes; byside < (TGeo::NSides); byside++)
942  {
943  if (refel->SideIsUndefined(byside)) {
944  continue;
945  }
946  TPZGeoElSide ElemSide(refel,byside);
947  TPZGeoElSide NextSide(ElemSide.Neighbour());
948  if(!NextSide.Element()) continue;
949  while(NextSide.Element() != ElemSide.Element())
950  {
951  if(NextSide.Exists() && !NextSide.Element()->IsLinearMapping() && !NextSide.Element()->IsGeoBlendEl())
952  {
953  TPZGeoElSide NeighSide = NextSide;
954  TPZTransform<> NeighTransf(NeighSide.Dimension(),NeighSide.Dimension());
955  ElemSide.SideTransform3(NeighSide,NeighTransf);
956  SetNeighbourInfo(byside,NeighSide,NeighTransf);
957  break;
958  }
959  NextSide = NextSide.Neighbour();
960  }
961  }
962 }
963 
964 /*
965  //inseri este método pois o compilador não encontrou
966  template <class TGeo>
967  void TPZGeoBlend<TGeo>::Initialize(TPZVec<int> &nodeindexes)
968  {
969  TGeo::Initialize(nodeindexes);
970  }
971  */
972 
973 #include "TPZGeoCube.h"
974 #include "TPZGeoLinear.h"
975 #include "pzgeoquad.h"
976 #include "pzgeotriangle.h"
977 #include "pzgeoprism.h"
978 #include "pzgeotetrahedra.h"
979 #include "pzgeopyramid.h"
980 #include "pzgeopoint.h"
981 #include "pzcmesh.h"
982 
983 
984 using namespace pzgeom;
985 
986 // template <class TGeo>
987 // TPZGeoEl *pzgeom::TPZGeoBlend<TGeo>::CreateBCGeoEl(TPZGeoEl *orig, int side,int bc)
988 // {
989 // TPZStack<int> LowAllSides;
990 // TGeo::LowerDimensionSides(side,LowAllSides);
991 // LowAllSides.Push(side);
992 // if(side < 0 || side > TGeo::NSides-1)
993 // {
994 // DebugStop();
995 // return 0;
996 // }
997 // bool straight = true;
998 // TPZGeoMesh *gmesh = orig->Mesh();
999 // for(int lowside = 0; lowside < LowAllSides.NElements(); lowside++)
1000 // {
1001 // if(LowAllSides[lowside] >= TGeo::NNodes && Neighbour(LowAllSides[lowside],gmesh).Element())
1002 // {
1003 // straight = false;
1004 // }
1005 // }
1006 // if(straight)
1007 // {
1008 // return TGeo::CreateBCGeoEl(orig,side,bc);
1009 // }
1010 // else
1011 // {
1012 // TPZGeoEl *newel = CreateBCGeoBlendEl(orig,side,bc);
1013 // return newel;
1014 // }
1015 
1016 // }
1017 
1018 // template <class TGeo>
1019 // TPZGeoEl *pzgeom::TPZGeoBlend<TGeo>::CreateBCGeoBlendEl(TPZGeoEl *orig,int side,int bc)
1020 // {
1021 // int ns = orig->NSideNodes(side);
1022 // TPZManVector<int64_t> nodeindices(ns);
1023 // int in;
1024 // for(in=0; in<ns; in++)
1025 // {
1026 // nodeindices[in] = orig->SideNodeIndex(side,in);
1027 // }
1028 // int64_t index;
1029 
1030 // TPZGeoMesh *mesh = orig->Mesh();
1031 // MElementType type = orig->Type(side);
1032 
1033 // TPZGeoEl *newel = mesh->CreateGeoBlendElement(type, nodeindices, bc, index);
1034 // TPZGeoElSide me(orig,side);
1035 // TPZGeoElSide newelside(newel,newel->NSides()-1);
1036 
1037 // newelside.InsertConnectivity(me);
1038 // newel->Initialize();
1039 
1040 // return newel;
1041 // }
1042 
1043 
1047 // template <class TGeo>
1048 // TPZGeoEl *pzgeom::TPZGeoBlend<TGeo>::CreateGeoElement(TPZGeoMesh &mesh, MElementType type,
1049 // TPZVec<int64_t>& nodeindexes,
1050 // int matid,
1051 // int64_t& index)
1052 // {
1053 // return CreateGeoElementMapped(mesh,type,nodeindexes,matid,index);
1054 // }
1055 
1057 /* @param gmesh mesh in which the element should be inserted
1058  @param matid material id of the element
1059  @param lowercorner (in/out) on input lower corner o the cube where the element should be created, on exit position of the next cube
1060  @param size (in) size of space where the element should be created
1061  */
1062 
1063 #include "TPZWavyLine.h"
1064 
1065 template <class TGeo>
1066 void pzgeom::TPZGeoBlend<TGeo>::InsertExampleElement(TPZGeoMesh &gmesh, int matid, TPZVec<REAL> &lowercorner, TPZVec<REAL> &size)
1067 {
1068 
1069  TGeo::InsertExampleElement(gmesh, -1, lowercorner, size);
1070  int64_t elid = gmesh.ElementVec().NElements()-1;
1071  TPZManVector<int64_t,3> nodeindexes(8);
1072  TPZGeoEl *gel = gmesh.Element(elid);
1073  int NNodes = TGeo::NCornerNodes;
1074  for (int i=0; i<NNodes; i++) {
1075  nodeindexes[i] = gel->NodeIndex(i);
1076  }
1077  TPZGeoElRefPattern<TPZWavyLine> *gelwave = new TPZGeoElRefPattern<TPZWavyLine>(nodeindexes,matid,gmesh);
1078  TPZManVector<REAL,3> wavedir(3,0.02);
1079  wavedir[0] = 0.;
1080  gelwave->Geom().SetData(wavedir, 2);
1081  delete gel;
1082  int64_t index;
1083  gmesh.CreateGeoBlendElement(TGeo::Type(), nodeindexes, matid, index);
1084 }
1085 
1086 
1087 
1088 
1089 template class pzgeom::TPZGeoBlend<TPZGeoCube>;
1090 template class pzgeom::TPZGeoBlend<TPZGeoTriangle>;
1091 template class pzgeom::TPZGeoBlend<TPZGeoPrism>;
1092 template class pzgeom::TPZGeoBlend<TPZGeoPyramid>;
1093 template class pzgeom::TPZGeoBlend<TPZGeoTetrahedra>;
1094 template class pzgeom::TPZGeoBlend<TPZGeoQuad>;
1095 template class pzgeom::TPZGeoBlend<TPZGeoLinear>;
1096 template class pzgeom::TPZGeoBlend<TPZGeoPoint>;
1097 
1098 #ifdef _AUTODIFF
1099 
1101 #define IMPLEMENTBLEND(TGEO,CREATEFUNCTION) \
1102 \
1103 template void pzgeom::TPZGeoBlend<TGEO>::X(TPZFMatrix<REAL> &coord, TPZVec<REAL> &par, TPZVec<REAL> &result) const; \
1104 template void pzgeom::TPZGeoBlend<TGEO>::X(TPZFMatrix<REAL> &coord, TPZVec<Fad<REAL>> &par, TPZVec<Fad<REAL>> &result) const; \
1105 template void pzgeom::TPZGeoBlend<TGEO>::GradX(TPZFMatrix<REAL> &coord, TPZVec<REAL> &xiInterior, TPZFMatrix<REAL> &gradx) const;\
1106 template void pzgeom::TPZGeoBlend<TGEO>::GradX(TPZFMatrix<REAL> &coord, TPZVec<Fad<REAL>> &xiInterior, TPZFMatrix<Fad<REAL>> &gradx) const;\
1107 template class \
1108 TPZRestoreClass< TPZGeoElRefPattern<TPZGeoBlend<TGEO> >>; \
1109 \
1110 template class TPZGeoElRefLess<TPZGeoBlend<TGEO> >;\
1111 template class TPZGeoElRefPattern<TPZGeoBlend<TGEO> >;
1112 
1113 IMPLEMENTBLEND(pzgeom::TPZGeoPoint,CreatePointEl)
1114 IMPLEMENTBLEND(pzgeom::TPZGeoLinear,CreateLinearEl)
1115 IMPLEMENTBLEND(pzgeom::TPZGeoQuad,CreateQuadEl)
1116 IMPLEMENTBLEND(pzgeom::TPZGeoTriangle,CreateTriangleEl)
1117 IMPLEMENTBLEND(pzgeom::TPZGeoCube,CreateCubeEl)
1118 IMPLEMENTBLEND(pzgeom::TPZGeoPrism,CreatePrismEl)
1119 IMPLEMENTBLEND(pzgeom::TPZGeoPyramid,CreatePyramEl)
1120 IMPLEMENTBLEND(pzgeom::TPZGeoTetrahedra,CreateTetraEl)
1121 
1122 #undef IMPLEMENTBLEND
1123 
1124 #else
1125 
1127 #define IMPLEMENTBLEND(TGEO,CREATEFUNCTION) \
1128 \
1129 template void pzgeom::TPZGeoBlend<TGEO>::X(TPZFMatrix<REAL> &coord, TPZVec<REAL> &par, TPZVec<REAL> &result) const; \
1130 template void pzgeom::TPZGeoBlend<TGEO>::GradX(TPZFMatrix<REAL> &coord, TPZVec<REAL> &xiInterior, TPZFMatrix<REAL> &gradx) const;\
1131 template class \
1132 TPZRestoreClass< TPZGeoElRefPattern<TPZGeoBlend<TGEO> >>; \
1133 \
1134 template class TPZGeoElRefLess<TPZGeoBlend<TGEO> >;\
1135 template class TPZGeoElRefPattern<TPZGeoBlend<TGEO> >;
1136 
1137 IMPLEMENTBLEND(pzgeom::TPZGeoPoint,CreatePointEl)
1138 IMPLEMENTBLEND(pzgeom::TPZGeoLinear,CreateLinearEl)
1139 IMPLEMENTBLEND(pzgeom::TPZGeoQuad,CreateQuadEl)
1140 IMPLEMENTBLEND(pzgeom::TPZGeoTriangle,CreateTriangleEl)
1141 IMPLEMENTBLEND(pzgeom::TPZGeoCube,CreateCubeEl)
1142 IMPLEMENTBLEND(pzgeom::TPZGeoPrism,CreatePrismEl)
1143 IMPLEMENTBLEND(pzgeom::TPZGeoPyramid,CreatePyramEl)
1144 IMPLEMENTBLEND(pzgeom::TPZGeoTetrahedra,CreateTetraEl)
1145 
1146 #undef IMPLEMENTBLEND
1147 
1148 
1149 #endif
1150 
1151 
1152 
1153 /*@orlandini : I REALLY dont know why is this here, so I have commented the following lines.
1154 If it breaks something, I am sorry.*/
1155 
1156 //#include "pznoderep.h.h"
1157 //template class pzgeom::TPZNodeRep<8,TPZGeoBlend<TPZGeoCube> >;
1158 //template class pzgeom::TPZNodeRep<6,TPZGeoBlend<TPZGeoPrism> >;
1159 
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
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
pzgeoelside.h
Contains declaration of TPZGeoElSide class which represents an element and its side, and TPZGeoElSideIndex class which represents an TPZGeoElSide index.
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.
pzgeoelrefless.h
Contains declaration of TPZGeoElRefLess class which implements the mapping between the master element...
TPZWavyLine.h
CreateQuadEl
TPZCompEl * CreateQuadEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational quadrilateral element.
Definition: pzcreateapproxspace.cpp:88
TPZFMatrix::GramSchmidt
void GramSchmidt(TPZFMatrix< TVar > &Orthog, TPZFMatrix< TVar > &TransfToOrthog)
This method implements a Gram Schimidt method. this = Orthog.TransfToOrthog.
Definition: pzfmatrix.cpp:341
pzgeom::TPZGeoBlend
Implements a blending map from curved boundaries to the interior of the element. Geometry.
Definition: tpzgeoblend.h:27
CreatePrismEl
TPZCompEl * CreatePrismEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational prismal element.
Definition: pzcreateapproxspace.cpp:109
MElementType_Name
std::string MElementType_Name(MElementType elType)
Returns the name of the element type.
Definition: pzeltype.h:127
TPZGeoElSide
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZChunkVector::NElements
int64_t NElements() const
Access method to query the number of elements of the vector.
Definition: TPZChunkVector.h:236
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
LOGPZ_WARN
#define LOGPZ_WARN(A, B)
Define log for warnings.
Definition: pzlog.h:91
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
pzgeom::TPZGeoLinear
Implements the geometry of a one dimensional linear element. Geometry.
Definition: TPZGeoLinear.h:30
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::SideTransform3
void SideTransform3(TPZGeoElSide neighbour, TPZTransform<> &t)
Accumulates the transformations from the current element/side to the neighbour/side.
Definition: pzgeoelside.cpp:601
TPZGeoMesh::CreateGeoBlendElement
virtual TPZGeoEl * CreateGeoBlendElement(MElementType type, TPZVec< int64_t > &nodeindexes, int matid, int64_t &index)
Creates a geometric element in same fashion of CreateGeoElement but here the elements are blend...
Definition: pzgmesh.cpp:1400
pzgeom::TPZGeoBlend::X
void X(TPZFMatrix< REAL > &cornerco, TPZVec< T > &par, TPZVec< T > &result) const
Get the coordinates of the point at geometric elements from coordinates of the parametric point at th...
Definition: tpzgeoblend.cpp:486
CreatePointEl
TPZCompEl * CreatePointEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational point element.
Definition: pzcreateapproxspace.cpp:71
pzgeom::TPZGeoBlend::Print
void Print(std::ostream &out=std::cout) const
Print all relevant data of the element to cout.
Definition: tpzgeoblend.cpp:925
tpzgeoelmapped.h
Contains declaration of TPZGeoElMapped class which implements a geometric element using its ancestral...
pzgeom::TPZGeoBlend::Initialize
void Initialize(TPZGeoEl *refel)
Initialize the element checking connectivity on all sides.
Definition: tpzgeoblend.cpp:938
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
CreateLinearEl
TPZCompEl * CreateLinearEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational linear element.
Definition: pzcreateapproxspace.cpp:79
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
TPZGeoLinear.h
TPZGeoElRefPattern
Implements a generic geometric element which is refined according to a generic refinement pattern...
Definition: pzgmesh.h:35
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
tpzgeoblend.h
Contains the TPZGeoBlend class which implements a blending map from curved boundaries to the interior...
TPZGeoMesh::Element
TPZGeoEl * Element(int64_t iel)
Definition: pzgmesh.h:139
pzgeotetrahedra.h
Contains the TPZGeoTetrahedra class which implements the geometry of a tetrahedral element...
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
MElementType_NNodes
int MElementType_NNodes(MElementType elType)
constant which defines the type of HDiv approximation space
Definition: pzeltype.h:80
pzgeom::TPZGeoBlend::ResetBlendConnectivity
bool ResetBlendConnectivity(const int64_t &side, const int64_t &index)
Definition: tpzgeoblend.cpp:48
pzgeom::TPZWavyLine::SetData
void SetData(TPZVec< REAL > &wavedir, int numwaves)
Definition: TPZWavyLine.h:58
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
Shape
void Shape(TPZVec< REAL > &pt, TPZVec< int > orders, TPZVec< TPZTransform< REAL > > &transvec, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
Definition: pzgenericshape.h:32
pzgeom::TPZGeoBlend::InsertExampleElement
static void InsertExampleElement(TPZGeoMesh &gmesh, int matid, TPZVec< REAL > &lowercorner, TPZVec< REAL > &size)
Method which creates a geometric boundary condition element based on the current geometric element...
Definition: tpzgeoblend.cpp:1066
pzgeom::TPZGeoQuad
Implements the geometry of a quadrilateral element. Geometry.
Definition: pzgeoquad.h:26
CreateTetraEl
TPZCompEl * CreateTetraEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational tetrahedral element.
Definition: pzcreateapproxspace.cpp:121
tpzgeoelrefpattern.h
Contains declaration of TPZGeoElRefPattern class which implements a generic geometric element which i...
pzgeom::TPZGeoBlend::GradX
void GradX(TPZFMatrix< REAL > &cornerco, TPZVec< T > &par, TPZFMatrix< T > &gradx) const
Computes the gradient of the transformation for parametric point at master element.
Definition: tpzgeoblend.cpp:84
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZExtractVal::val
static REAL val(const int number)
Definition: pzextractval.h:21
pzcmesh.h
Contains declaration of TPZCompMesh class which is a repository for computational elements...
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
pzgeom::TPZGeoPyramid
Implements the geometry of pyramid element. Geometry.
Definition: pzgeopyramid.h:24
CreatePyramEl
TPZCompEl * CreatePyramEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational pyramidal element.
Definition: pzcreateapproxspace.cpp:115
TPZGeoCube.h
Contains the TPZGeoCube class which implements the geometry of hexahedra element. ...
pzgeopoint.h
Contains the TPZGeoPoint class which implements the geometry of a point element or 0-D element...
pzgeom::TPZGeoPoint
Implements the geometry of a point element. Geometry.
Definition: pzgeopoint.h:39
pzgeom::TPZGeoTriangle
Implements the geometry of a triangle element. Geometry.
Definition: pzgeotriangle.h:28
pzgeom::TPZGeoBlend::Jacobian
void Jacobian(TPZFMatrix< REAL > &cornerco, TPZVec< REAL > &par, TPZFMatrix< REAL > &jacobian, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Computes the Jacobian for parametric point at master element.
Definition: tpzgeoblend.cpp:733
TPZGeoElSide::Element
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
CreateCubeEl
TPZCompEl * CreateCubeEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational cube element.
Definition: pzcreateapproxspace.cpp:103
pzgeom::TPZGeoTetrahedra
Implements the geometry of a tetrahedral element. Geometry.
Definition: pzgeotetrahedra.h:26
TPZGeoEl::SideIsUndefined
virtual int SideIsUndefined(int side)=0
Returns 1 if the side has not been defined by buildconnectivity.
TPZGeoEl::Id
int64_t Id() const
Returns the Id of the element.
Definition: pzgeoel.h:223
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
MElementType
MElementType
Define the element types.
Definition: pzeltype.h:52
TPZGeoElSide::Dimension
int Dimension() const
the dimension associated with the element/side
Definition: pzgeoelside.cpp:593
pzgeom::TPZGeoBlend::SetNeighbourInfo
void SetNeighbourInfo(int side, TPZGeoElSide &neigh, TPZTransform<> &trans) override
Definition: tpzgeoblend.cpp:57
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
pzgeom::TPZGeoPrism
Implements the geometry of a prism element. Geometry.
Definition: pzgeoprism.h:27
TPZVec::Fill
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
TPZGeoElSide::Exists
int Exists() const
Definition: pzgeoelside.h:175
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
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.
pzgeom::TPZGeoBlend::IsLinearMapping
bool IsLinearMapping(int side) const
Definition: tpzgeoblend.cpp:20
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
IMPLEMENTBLEND
#define IMPLEMENTBLEND(TGEO, CREATEFUNCTION)
CreateGeoElement -> TPZGeoBlend.
Definition: tpzgeoblend.cpp:1127
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
pzgeoprism.h
Contains the TPZGeoPrism class which implements the geometry of a prism element.
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
ZeroTolerance
REAL ZeroTolerance()
Returns the tolerance to Zero value. Actually: .
Definition: pzreal.h:633
pzgeom::TPZGeoCube
Implements the geometry of hexahedra element. Geometry.
Definition: TPZGeoCube.h:28
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
CreateTriangleEl
TPZCompEl * CreateTriangleEl(TPZGeoEl *gel, TPZCompMesh &mesh, int64_t &index)
Creates computational triangular element.
Definition: pzcreateapproxspace.cpp:97
TPZGeoMesh::ElementVec
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138
Generated on Sun Aug 16 2020 11:17:01 for NeoPZ by   doxygen 1.8.13