NeoPZ
TPZRefPatternTools.cpp
Go to the documentation of this file.
1 
5 /*
6  * Created by Cesar Lucci on 10/03/10.
7  * Copyright 2010 LabMeC. All rights reserved.
8  *
9  */
10 
11 #include "TPZRefPatternTools.h"
12 #include "TPZVTKGeoMesh.h"
13 #include "TPZRefPattern.h"
14 #include "TPZRefPatternDataBase.h"
15 
16 using namespace std;
17 
19 {
20 
21 }
22 
24 {
25 
26 }
27 
29  std::list<TPZAutoPointer<TPZRefPattern> > &refs)
30 {
31  if(!gel) return;
32 
33  refs.clear();
34  if(gel->HasSubElement())
35  {
36  refs.push_back(gel->GetRefPattern());
37  return;
38  }
39 
40  // first we build the refinement patterns associated with the neighbours of the current element
41  int side, nsides, nnodes;
42  nsides = gel->NSides();
43  nnodes = gel->NCornerNodes();
44  TPZManVector<TPZAutoPointer<TPZRefPattern>, 27> NeighSideRefPatternVec(nsides,0);
45 
46  for(side = nnodes; side < nsides; side++)
47  {
48  TPZGeoElSide gelside(gel, side);
49  TPZGeoElSide neighside = gelside.Neighbour();
50  if(!neighside.Element()){
51  DebugStop();
52  }
53  while(neighside != gelside)
54  {
55  if(neighside.NSubElements() > 1)
56  {
57  TPZAutoPointer<TPZRefPattern> neighRefp = neighside.Element()->GetRefPattern();
58  if(neighRefp)
59  {
60  TPZTransform<> trans = neighside.NeighbourSideTransform(gelside);
61  NeighSideRefPatternVec[side] = neighRefp->SideRefPattern(neighside.Side(),trans);
62  break;
63  }
64  }
65  neighside = neighside.Neighbour();
66  }
67  }
68 
69  // having the refinement patterns associated with the sides, look for compatible refinement patterns
70  std::list< TPZAutoPointer<TPZRefPattern> > gelReflist = gRefDBase.RefPatternList(gel->Type());
71  std::list< TPZAutoPointer<TPZRefPattern> >::iterator gelReflistIt;
72 
73  for(gelReflistIt = gelReflist.begin(); gelReflistIt != gelReflist.end(); gelReflistIt++)
74  {
75  // compare the side refinement patterns
76  for(side = nnodes; side < nsides; side++)
77  {
78  TPZAutoPointer<TPZRefPattern> GelSideRefPattern = (*gelReflistIt)->SideRefPattern(side);
79  TPZAutoPointer<TPZRefPattern> NeighSideRefPattern = NeighSideRefPatternVec[side];
80 
81  if(GelSideRefPattern && NeighSideRefPattern)
82  {
83  if( !( (TPZRefPattern &)GelSideRefPattern == NeighSideRefPattern) )
84  {
85  break;
86  }
87  }
88  }
89 
90  // if all refinement patterns are equal
91  if(side == nsides)
92  {
93  refs.push_back(*gelReflistIt);
94  }
95  }
96 }
97 
99  std::map<int, std::pair<TPZGeoEl *,
100  std::map<int,int> > > &neighCorresp)
101 {
102  neighCorresp.clear();
103  TPZAutoPointer<TPZRefPattern> modelPat;//NULL
104 
105  int nnodes = gel->NNodes();
106  int nsides = gel->NSides();
107 
108  bool neigh_isrefined, topolAreCompatibles;
109 
110  std::list< TPZAutoPointer<TPZRefPattern> > candidates = gRefDBase.RefPatternList(gel->Type());
111  std::list< TPZAutoPointer<TPZRefPattern> >::iterator it;
112 
113  int side;
114  for(it = candidates.begin(); it != candidates.end(); it++)
115  {
116  modelPat = *it;
117  for(side = nnodes; side < nsides; side++)
118  {
119  TPZGeoElSide gelside(gel, side);
120  TPZGeoElSide neighside = gelside.Neighbour();
121 
122  neigh_isrefined = false;
123  topolAreCompatibles = false;
124 
125  //Se meu lado eh aresta e nao tem vizinho, significa que nao precisa ser refinado!
126  //Mas se o candidato a modelRefPattern quer refinar por este lado, jah nao serve... vamos para outro candidato!
127  if(neighside == gelside && gelside.Dimension() == 1 && modelPat->SideRefPattern(side) && side != nsides-1)
128  {
129  neigh_isrefined = true;
130  topolAreCompatibles = false;
131  break;
132  }
133  //Se eu nao tenho vizinho pelo lado (nsides-1), o candidato a modelRefPattern pode refina-lo!!!
134  else if(neighside == gelside && side == nsides-1)
135  {
136  continue;
137  }
138 
139  while(neighside != gelside)
140  {
141  if(neighside.NSubElements() > 1)
142  {
143  neigh_isrefined = true;
144  TPZAutoPointer<TPZRefPattern> neighRefp = neighside.Element()->GetRefPattern()->SideRefPattern(neighside.Side());
145  TPZAutoPointer<TPZRefPattern> modelsideRefp = modelPat->SideRefPattern(side);
146 
147  if(neighRefp && modelsideRefp)
148  {
149  TPZTransform<> transBetweenNeigh = neighside.NeighbourSideTransform(gelside);
150  TPZTransform<> transBetweenSides = gel->SideToSideTransform(side, nsides-1);
151 
152  std::map<int, int> pairedNodes;
153 
154  //correspondencias entre nohs entre LADO DO PADRAO MODELO e LADO DO PADRAO VIZINHO
155  topolAreCompatibles = CompareTopologies(neighRefp, modelsideRefp, transBetweenNeigh, pairedNodes);
156 
157  if(topolAreCompatibles)
158  {
159  neighCorresp[side] = make_pair(neighside.Element(), pairedNodes);
160  break;
161  }
162  }
163  }
164  neighside = neighside.Neighbour();
165  }
166  if( (neigh_isrefined && !topolAreCompatibles) || (gelside.Dimension() == 1 && !neigh_isrefined && modelPat->SideRefPattern(side)) )
167  {
168  neighCorresp.clear();
169  break;
170  }
171  }
172  if(side == nsides)
173  {
174  return modelPat;
175  }
176  }
177 
178  return NULL;
179 }
180 
182  TPZVec<int> &sidestorefine)
183 {
184  if(!gel)
185  {
186  return NULL;
187  }
188 
189  std::list<TPZAutoPointer<TPZRefPattern> > patlist, perfectmatch;
191 
192  // totototo
193  //std::cout << "gel index " << gel->Index() << " sides " << sidestorefine << std::endl;
194  std::list<TPZAutoPointer<TPZRefPattern> >::iterator it;
195  for(it = patlist.begin(); it != patlist.end(); it++)
196  {
197  if( !(*it) )
198  {
199  continue;
200  }
201  TPZGeoEl * gelTemp = (*it)->Element(0);
202  int ncorners = gelTemp->NCornerNodes();
203  int nsides = gelTemp->NSides();
204  if(gelTemp->Dimension() == 3)
205  {
206  nsides--;
207  }
208 
209  int is;
210  for(is = ncorners; is < nsides; is++)
211  {
212  if( sidestorefine[is] != (*it)->NSideNodes(is) )
213  {
214  break;
215  }
216  }
217  if(is == nsides)
218  {
219  perfectmatch.push_back(*it);
220  }
221  }
222  //std::cout << "perfect match size " << perfectmatch.size() << std::endl;
223 
224  patlist.clear();
225  if (perfectmatch.size() == 1) {
226  return *perfectmatch.begin();
227  } else if(perfectmatch.size() > 1)
228  {
229  // return the refpattern with the least number of elements
230 
231  int minsubel = 30;
232  for (auto it=perfectmatch.begin(); it!= perfectmatch.end(); it++) {
233  TPZAutoPointer<TPZRefPattern> refpat = *it;
234  if (refpat->NSubElements() < minsubel) {
235  minsubel = refpat->NSubElements();
236  patlist.clear();
237  patlist.push_back(*it);
238  }
239  else if(refpat->NSubElements() == minsubel)
240  {
241  patlist.push_back(*it);
242  }
243  }
244  if (patlist.size() != 1) {
245  for (auto it2=patlist.begin(); it2 != patlist.end(); it2++) {
246  (*it2)->fRefPatternMesh.Print();
247  }
248  DebugStop();
249  }
250  return *patlist.begin();
251  }
252  else
253  {
254  return NULL;
255  }
256 }
257 
259 {
260  //Criando uma malha local a partir do vetor de elementos.
261  TPZGeoMesh refGMesh;
262  GenerateGMeshFromElementVec(realMeshElementVec,refGMesh);
263 
264  TPZAutoPointer<TPZRefPattern> refp = new TPZRefPattern(refGMesh);
266  if(refpFound)
267  {
268  ModifyElementsBasedOnRefpFound(refpFound,refp,realMeshElementVec);
269  return refpFound;
270  }
271  else
272  {
274  return refp;
275  }
276 
277  return NULL;
278 }
279 
281 {
282  std::map<int64_t,int64_t> gl2lcNdIdx;//global to local node index
283  std::map<int64_t,int64_t> gl2lcElIdx;//global to local element index
284 
285  refGMesh.ElementVec().Resize(elementVec.NElements());
286  for(int el = 0; el < elementVec.NElements(); el++)
287  {
288  TPZGeoEl *gel = elementVec[el];
289  gl2lcElIdx[gel->Index()] = el;
290  }
291  for(int el = 0; el < elementVec.NElements(); el++)
292  {
293  TPZGeoEl *gel = elementVec[el];
294  if(!gel) DebugStop();
295  for(int n = 0; n < gel->NCornerNodes(); n++)
296  {
297  int nodeindex = gel->NodeIndex(n);
298  if(gl2lcNdIdx.find(nodeindex) != gl2lcNdIdx.end())
299  {
300  continue;
301  }
302 
303  TPZManVector<REAL,3> sonNodeCoord(3);
304  gel->NodePtr(n)->GetCoordinates(sonNodeCoord);
305 
306  int oldSize = refGMesh.NNodes();
307  refGMesh.NodeVec().Resize(oldSize+1);
308  refGMesh.NodeVec()[oldSize].SetCoord(sonNodeCoord);
309  refGMesh.NodeVec()[oldSize].SetNodeId(oldSize);
310 
311  gl2lcNdIdx[nodeindex] = oldSize;
312  }
313  gel->ClonePatchEl(refGMesh,gl2lcNdIdx,gl2lcElIdx);
314  }
315 }
316 
319  TPZVec<TPZGeoEl *> &elementVec)
320 {
321  TPZGeoMesh & refpGMesh = refp->RefPatternMesh();
322  TPZGeoMesh & refpFoundGMesh = refpFound->RefPatternMesh();
323 
324  //pareando nohs
325  std::map<int64_t,int64_t> Node_RefpFound2Refp;
326  for(int64_t n = 0; n < refpGMesh.NNodes(); n++)
327  {
328  TPZManVector<REAL,3> coord(3);
329  refpGMesh.NodeVec()[n].GetCoordinates(coord);
330 
331  int nodeIndex;
332  TPZGeoNode * node = refpFoundGMesh.FindNode(coord,nodeIndex);
333  if(!node)
334  {
335  DebugStop();
336  }
337 
338  Node_RefpFound2Refp[nodeIndex] = n;
339  }
340 
341  TPZManVector<TPZGeoEl *> subels(elementVec.size()-1,0);
342  for(int64_t el=1; el < refpGMesh.NElements(); el++)
343  {
344  TPZManVector<int64_t> nodeIndicesVecRefp;
345  std::set<int64_t> nodeindicesRefp;
346  refpGMesh.ElementVec()[el]->GetNodeIndices(nodeIndicesVecRefp);
347  if(nodeIndicesVecRefp.NElements())
348  {
349  int sz = nodeIndicesVecRefp.size();
350  int64_t *zero = &nodeIndicesVecRefp[0];
351  nodeindicesRefp.insert(zero,zero+sz);
352  }
353  else
354  {
355  DebugStop();
356  }
357 
358  //pareando subelemento
359  int64_t el2 = 1;
360  TPZManVector<int64_t> nodeIndicesVecRefpFound;
361  for(el2 = 1; el2 < refpFoundGMesh.NElements(); el2++)
362  {
363  refpFoundGMesh.ElementVec()[el2]->GetNodeIndices(nodeIndicesVecRefpFound);
364  std::set<int64_t> nodeindicesRefpFound;
365  for(int n = 0; n < nodeIndicesVecRefpFound.NElements(); n++)
366  {
367 #ifdef PZDEBUG
368  if(Node_RefpFound2Refp.find(nodeIndicesVecRefpFound[n]) == Node_RefpFound2Refp.end())
369  {
370  DebugStop();
371  }
372 #endif
373 
374  nodeindicesRefpFound.insert(Node_RefpFound2Refp[nodeIndicesVecRefpFound[n]]);
375  }
376  if(nodeindicesRefp == nodeindicesRefpFound)
377  {
378  break;
379  }
380  }
381  if(el2 == refpFoundGMesh.NElements())
382  {
383  DebugStop();
384  }
385  subels[el2-1] = elementVec[el];
386 
387  //Elementos pareados
388  TPZManVector<int64_t,8> newTopol(refpFoundGMesh.ElementVec()[el2]->NNodes());
389  for(int64_t n = 0; n < refpFoundGMesh.ElementVec()[el2]->NNodes(); n++)
390  {
391  std::map<int64_t,int64_t>::iterator it = Node_RefpFound2Refp.find(refpFoundGMesh.ElementVec()[el2]->NodeIndex(n));
392  if(it == Node_RefpFound2Refp.end())
393  {
394  DebugStop();
395  }
396  int newNodeIndex = -1;
397  for(int nn = 0; nn < refpGMesh.ElementVec()[el]->NNodes(); nn++)
398  {
399  if(refpGMesh.ElementVec()[el]->NodeIndex(nn) == it->second)
400  {
401  newNodeIndex = nn;
402  break;
403  }
404  }
405  if(newNodeIndex == -1)
406  {
407  DebugStop();
408  }
409  newTopol[n] = elementVec[el]->NodeIndex(newNodeIndex);
410  }
411  for(int64_t n = 0; n < elementVec[el]->NNodes(); n++)
412  {
413  elementVec[el]->SetNodeIndex(n,newTopol[n]);
414  }
415  }
416 
417  for(int64_t el = 1; el < elementVec.NElements(); el++)
418  {
419  elementVec[el] = subels[el-1];
420  }
421  //Lembre-se de chamar elementVec[0]->Mesh()->BuildConnectivity(), uma vez que houve alteracao
422  //das declaracoes topologicas dos elementos do elementVec, modificando conectividades (vizinhanca).
423 }
424 
426 {
427  if(!gel)
428  {
429  return NULL;
430  }
431  TPZManVector<int,27> sidestorefine;
432 
433  bool thereIsAnyNeighbourRefined = TPZRefPatternTools::SidesToRefine(gel, sidestorefine);
434  if(!thereIsAnyNeighbourRefined)
435  {
436  return 0;
437  }
438  TPZAutoPointer<TPZRefPattern> pat = PerfectMatchRefPattern(gel, sidestorefine);
439 
440  if(!pat)
441  {
442  //if no refpattern matches perfectly...
443  std::map<int, std::pair<TPZGeoEl *, std::map<int,int> > > neighCorresp;
445  if(modelPat)
446  {
447  //modelPat->Print();
448  pat = TPZRefPatternTools::DragModelPatNodes(gel, modelPat, neighCorresp);
449  }
450  else
451  {
452  std::cout << "\nModelRefPattern NOT FOUND in " << __PRETTY_FUNCTION__ << std::endl;
453  std::cout << "You should create it and add in Refinement Patterns Folder!" << std::endl;
454  std::cout << "Open file ModelRefPatternNOTFOUND.vtk in Paraview to see the neighbourhood\n";
455 
456  std::ofstream outNotFound("ModelRefPatternNOTFOUND.vtk");
457  TPZVTKGeoMesh::PrintGMeshVTKneighbourhood(gel->Mesh(), gel->Index(), outNotFound);
458 
459  DebugStop();
460  }
461  }
462 
463  return pat;
464 }
465 
468  std::map<int, std::pair<TPZGeoEl *,
469  std::map<int,int> > > &neighCorresp)
470 {
471  TPZAutoPointer<TPZRefPattern> modelPat_copy = new TPZRefPattern(modelPat);
472  int nnodes = gel->NNodes();
473  int nsides = gel->NSides();
474 
475  for(int side = nnodes; side < nsides; side++)
476  {
477  std::map<int, std::pair<TPZGeoEl *, std::map<int,int> > >::iterator correspIt = neighCorresp.find(side);
478  if(correspIt == neighCorresp.end())
479  {
480  continue;
481  }
482  TPZGeoEl * neighStored = correspIt->second.first;
483  TPZGeoElSide gelside(gel, side);
484  TPZGeoElSide neighside = gelside.Neighbour();
485 
486  while(neighside != gelside)
487  {
488  if(neighside.Element() == neighStored)
489  {
490  //int SideOfNeighSide = neighside.Side();
491  TPZAutoPointer<TPZRefPattern> neighRefp = neighside.Element()->GetRefPattern()->SideRefPattern(neighside.Side());
492  TPZAutoPointer<TPZRefPattern> modelsideRefp = modelPat_copy->SideRefPattern(side);
493 
494  TPZTransform<> transBetweenNeigh = neighside.NeighbourSideTransform(gelside);
495  TPZTransform<> transBetweenSides = gel->SideToSideTransform(side, nsides-1);
496 
497  //correspondencias entre nohs entre LADO DO PADRAO MODELO e LADO DO PADRAO VIZINHO
498  std::map<int, int> pairedNodes = correspIt->second.second;
499 
500  //correspondencias entre nohs entre LADO DO PADRAO MODELO e o proprio PADRAO MODELO
501  std::map<int, int> pairSideNodes;
502  PairMeshesNodesMatchingCoordinates(modelsideRefp->RefPatternMesh(), modelPat_copy->RefPatternMesh(), transBetweenSides, pairSideNodes);
503  std::map<int, int>::iterator it, jt;
504  for(it = pairedNodes.begin(); it != pairedNodes.end(); it++)
505  {
506  TPZVec<REAL> NodeCoord(3,0.), sideNodeCoord(3,0.), neighNodeCoord(3,0.);
507 
508  neighRefp->RefPatternMesh().NodeVec()[it->first].GetCoordinates(neighNodeCoord);
509  transBetweenNeigh.Apply(neighNodeCoord, sideNodeCoord);
510  transBetweenSides.Apply(sideNodeCoord, NodeCoord);
511 
512  jt = pairSideNodes.find(it->second);
513  if(jt != pairSideNodes.end())
514  {
515  int gNode = jt->second;
516  for(int n = 0; n < 3; n++)
517  {
518  modelPat_copy->RefPatternMesh().NodeVec()[gNode].SetCoord(n, NodeCoord[n]);
519  }
520  }
521  }
522 
523  break;
524  }
525  neighside = neighside.Neighbour();
526  }
527  }
528  modelPat_copy->ComputeTransforms();
529  modelPat_copy->ComputePartition();
530  modelPat_copy->GenerateSideRefPatterns();
531  if(!gRefDBase.FindRefPattern(modelPat_copy))
532  {
533  gRefDBase.InsertRefPattern(modelPat_copy);
534  modelPat_copy->InsertPermuted();
535  }
536 
537  return modelPat_copy;
538 }
539 
542  TPZTransform<> &fromAtoB,
543  std::map<int, int> &pairedNodes)
544 {
545  pairedNodes.clear();
546 
547  if(refA->Type() != refB->Type())
548  {
549  //Their are NOT topologicaly compatibles
550  return false;
551  }
552 
553  TPZGeoMesh meshA = refA->RefPatternMesh();
554  TPZGeoMesh meshB = refB->RefPatternMesh();
555 
556  if(meshA.NNodes() != meshB.NNodes() || meshA.NElements() != meshB.NElements())
557  {
558  //Their are NOT topologicaly compatibles
559  return false;
560  }
561 
562  //Subgrouping sub-elements by types
563  std::map<MElementType, std::list<TPZGeoEl*> > eltypeA, eltypeB;
564  for(int el = 1; el < meshA.NElements(); el++)//Father is NOT included!
565  {
566  TPZGeoEl *gelA = meshA.ElementVec()[el];
567  TPZGeoEl *gelB = meshB.ElementVec()[el];
568 
569  eltypeA[gelA->Type()].push_back(gelA);
570  eltypeB[gelB->Type()].push_back(gelB);
571  }
572  if(eltypeA.size() != eltypeB.size())
573  {
574  //Their are NOT topologicaly compatibles
575  return false;
576  }
577  std::map<MElementType, std::list<TPZGeoEl*> >::iterator eltypeAit, eltypeBit;
578  for(eltypeAit = eltypeA.begin(); eltypeAit != eltypeA.end(); eltypeAit++)
579  {
580  eltypeBit = eltypeB.find(eltypeAit->first);
581  if( (eltypeBit == eltypeB.end()) || (eltypeAit->second.size() != eltypeBit->second.size()) )
582  {
583  //Their are NOT topologicaly compatibles
584  return false;
585  }
586  }
587 
588  //PARTY BEGINS!!!!!!!!!
589  PairMeshesCornerNodesMatchingCoordinates(meshA, meshB, fromAtoB, pairedNodes);
590 
591  TPZGeoEl* fatherA = meshA.ElementVec()[0];
592  int unpairedNNodes = meshA.NNodes() - pairedNodes.size();
593 
594 #ifdef PZDEBUG
595  if((int) pairedNodes.size() < fatherA->NNodes())
596  {
597  std::cout << "\nThere is something going wrong with meshA and/or meshB in " << __PRETTY_FUNCTION__ << std::endl;
598  std::cout << "father->ConnerNodes should be paired at least!!!\n\n";
599  DebugStop();
600  }
601 #endif
602 
603  std::map<int, int>::iterator pairedNodesIT;
604 
605  if(unpairedNNodes > 0)
606  {
607  for(int s = fatherA->NNodes(); s < fatherA->NSides(); s++)//pairing fatherA~fatherB (Edges)Nodes
608  {
609  TPZGeoElSide Aside(fatherA, s);
610  if(Aside.Dimension() == 1)
611  {
612  pairedNodesIT = pairedNodes.find(Aside.SideNodeIndex(0));
613  int nodeIniA = pairedNodesIT->first;
614  int nodeIniB = pairedNodesIT->second;
615 
616  pairedNodesIT = pairedNodes.find(Aside.SideNodeIndex(1));
617  int nodeFinA = pairedNodesIT->first;
618  int nodeFinB = pairedNodesIT->second;
619 
620  TPZManVector<int> NodesHuntedA;
621  NodesHunter(meshA, NodesHuntedA, nodeIniA, nodeFinA);
622 
623  TPZManVector<int> NodesHuntedB;
624  NodesHunter(meshB, NodesHuntedB, nodeIniB, nodeFinB);
625 
626  for(int n = 1; n < NodesHuntedA.size()-1; n++)
627  {
628  int nA = NodesHuntedA[n];
629  int nB = NodesHuntedB[n];
630 
631  if(pairedNodes.find(nA) == pairedNodes.end())
632  {
633  pairedNodes[nA] = nB;
634  unpairedNNodes--;
635  }
636  }
637  }
638  }
639  }
640 
641  //pairing father->(Faces&Volume)Nodes AND SUB-ELEMENTS!!!
642  int nRemainingSubels = refA->NSubElements();
643  TPZVec<bool> AsubelementIsAlreadyPaired(meshA.NElements(), false);
644  TPZVec<bool> BsubelementIsAlreadyPaired(meshB.NElements(), false);
645 
646 #define IwantPairedElementsToo
647 #ifdef IwantPairedElementsToo
648  std::map<int, int> pairedElements;
649  pairedElements[0] = 0;//fathers are already paired at this line of code
650 #endif
651 
652  std::list<TPZGeoEl*>::iterator elAit, elBit;
653  std::list< std::list<TPZGeoEl*>::iterator > pointedB;
654  TPZVec<int> Anodes(0);
655  TPZVec< TPZManVector<int,8> > Bnodes(0);
656  std::list< TPZVec<int> > BthatMatch;
657  int pairedNodesCount, unsuccessfulLoop = 0;
658  while(unsuccessfulLoop < 2 && nRemainingSubels > 0)
659  {
660  for(eltypeAit = eltypeA.begin(); eltypeAit != eltypeA.end(); eltypeAit++)
661  {
662  eltypeBit = eltypeB.find(eltypeAit->first);
663  for(elAit = eltypeAit->second.begin(); elAit != eltypeAit->second.end(); elAit++)
664  {
665  TPZGeoEl * A = *elAit;
666  Anodes.Resize(A->NNodes());
667  pairedNodesCount = 0;
668  for(int n = 0; n < A->NNodes(); n++)
669  {
670  Anodes[n] = A->NodeIndex(n);
671  pairedNodesIT = pairedNodes.find(Anodes[n]);
672  if(pairedNodesIT != pairedNodes.end())
673  {
674  pairedNodesCount++;
675  }
676  }
677  if(pairedNodesCount > 1)
678  {
679  BthatMatch.clear();
680  pointedB.clear();
681  for(elBit = eltypeBit->second.begin(); elBit != eltypeBit->second.end(); elBit++)
682  {
683  TPZGeoEl * B = *elBit;
685  for(int p = 0; p < Bnodes.NElements(); p++)
686  {
687  bool goodCandidate = true;
688  for(int n = 0; n < Anodes.NElements(); n++)
689  {
690  int nA = Anodes[n];
691  int nB = Bnodes[p][n];
692 
693  pairedNodesIT = pairedNodes.find(nA);
694  if(pairedNodesIT != pairedNodes.end() && pairedNodesIT->second != nB)
695  {
696  goodCandidate = false;
697  break;
698  }
699  }
700  if(goodCandidate)
701  {
702  BthatMatch.push_back(Bnodes[p]);
703  pointedB.push_back(elBit);
704  }
705  }
706  }//for elBit
707  }//if(pairedNodesCount > 1)
708  if(pairedNodesCount > 1 && BthatMatch.size() == 1)//se temos apenas 1 candidato da lista de B
709  {
710  for(int n = 0; n < Anodes.NElements(); n++)//pareando os nohs ainda nao pareados e deletando da lista os elementos pareados
711  {
712  pairedNodesIT = pairedNodes.find(Anodes[n]);
713  if(pairedNodesIT == pairedNodes.end() && unpairedNNodes > 0)
714  {
715  pairedNodes[Anodes[n]] = (*BthatMatch.begin())[n];
716  unpairedNNodes--;
717  }
718  }
719  nRemainingSubels--;
720 
721 #ifdef IwantPairedElementsToo
722  int Aid = (*(elAit))->Id();
723  int Bid = (*(*(pointedB.begin())))->Id();
724  pairedElements[Aid] = Bid;
725 #endif
726 
727  eltypeAit->second.erase(elAit);
728  eltypeBit->second.erase(*pointedB.begin());
729 
730  unsuccessfulLoop = 0;
731  break;
732  }
733  }//for elAit
734  }//for eltypeAit
735  unsuccessfulLoop++;
736  }//while
737 
738  if(nRemainingSubels == 0 && unpairedNNodes == 0)
739  {
740  //Their ARE topologicaly compatibles
741  return true;
742  }
743 
744  return false;
745 }
746 
748  TPZGeoMesh &meshB,
749  TPZTransform<> &fromAtoB,
750  std::map<int, int> &pairedNodes)
751 {
752  TPZVec<REAL> ANodeCoord(3,0.), BNodeCoord(3,0.);
753 
754  int Annodes = meshA.ElementVec()[0]->NNodes();
755  int Bnnodes = meshB.ElementVec()[0]->NNodes();
756  for(int nA = 0; nA < Annodes; nA++)
757  {
758  meshA.NodeVec()[meshA.ElementVec()[0]->NodeIndex(nA)].GetCoordinates(ANodeCoord);
759  fromAtoB.Apply(ANodeCoord, BNodeCoord);
760 
761  TPZVec<REAL> BNodeCoord_compare(3,0.);
762  for(int nB = 0; nB < Bnnodes; nB++)
763  {
764  double dif = 0.;
765  meshB.NodeVec()[meshB.ElementVec()[0]->NodeIndex(nB)].GetCoordinates(BNodeCoord_compare);
766  for(int c = 0; c < 3; c++)
767  {
768  dif += fabs(BNodeCoord[c] - BNodeCoord_compare[c]);
769  }
770  if(dif < 1.E-10)
771  {
772  int nA_Id = meshA.NodeVec()[meshA.ElementVec()[0]->NodeIndex(nA)].Id();
773  int nB_Id = meshB.NodeVec()[meshB.ElementVec()[0]->NodeIndex(nB)].Id();
774  pairedNodes[nA_Id] = nB_Id;
775 
776  break;
777  }
778  }
779  }
780 }
781 
783  TPZGeoMesh &meshB,
784  TPZTransform<> &fromAtoB,
785  std::map<int, int> &pairedNodes)
786 {
787  TPZVec<REAL> ANodeCoord(3,0.), BNodeCoord(3,0.);
788 
789  int Annodes = meshA.NNodes();
790  int Bnnodes = meshB.NNodes();
791  for(int nA = 0; nA < Annodes; nA++)
792  {
793  meshA.NodeVec()[nA].GetCoordinates(ANodeCoord);
794  fromAtoB.Apply(ANodeCoord, BNodeCoord);
795 
796  TPZVec<REAL> BNodeCoord_compare(3,0.);
797  for(int nB = 0; nB < Bnnodes; nB++)
798  {
799  double dif = 0.;
800  meshB.NodeVec()[nB].GetCoordinates(BNodeCoord_compare);
801  for(int c = 0; c < 3; c++)
802  {
803  dif += fabs(BNodeCoord[c] - BNodeCoord_compare[c]);
804  }
805  if(dif < 1.E-10)
806  {
807  int nA_Id = meshA.NodeVec()[nA].Id();
808  int nB_Id = meshB.NodeVec()[nB].Id();
809  pairedNodes[nA_Id] = nB_Id;
810 
811  break;
812  }
813  }
814  }
815 }
816 
818 {
819  std::string refpTypeName;
820  std::stringstream rpName;
821  TPZRefPattern *prefp = &refp;
822 
823  if(prefp == NULL)
824  {
825  std::cout << "Null refpattern parameter on " << __PRETTY_FUNCTION__ << std::endl;
826  return refpTypeName;
827  }
828 
829  TPZGeoEl *gel = refp.Element(0);
830  int ncorners = gel->NCornerNodes();
831  int nsides = gel->NSides();
832 
833  std::string perfix = gel->TypeName();
834  for(int i = 0; i < 3; i++)
835  {
836  rpName << perfix[i];
837  }
838  for(int s = 0; s < nsides; s++)
839  {
840  if(s < ncorners)
841  {
842  rpName << "0";
843  }
844  else
845  {
846  rpName << refp.NSideNodes(s);
847  }
848  }
849  rpName >> refpTypeName;
850 
851  if(refpTypeName.length() == 0)
852  {
853  refpTypeName = TPZRefPattern::fNonInitializedName;
854  }
855 
856  return refpTypeName;
857 }
858 
860 {
861  std::string refpTypeName;
862  std::stringstream rpName;
863 
864  if(!refp)
865  {
866  std::cout << "Null refpattern parameter on " << __PRETTY_FUNCTION__ << std::endl;
867  return refpTypeName;
868  }
869 
870  TPZGeoEl *gel = refp->Element(0);
871  int ncorners = gel->NCornerNodes();
872  int nsides = gel->NSides();
873 
874  std::string perfix = gel->TypeName();
875  for(int i = 0; i < 3; i++)
876  {
877  rpName << perfix[i];
878  }
879  for(int s = 0; s < nsides; s++)
880  {
881  if(s < ncorners)
882  {
883  rpName << "0";
884  }
885  else
886  {
887  rpName << refp->NSideNodes(s);
888  }
889  }
890  rpName >> refpTypeName;
891 
892  if(refpTypeName.length() == 0)
893  {
894  refpTypeName = TPZRefPattern::fNonInitializedName;
895  }
896 
897  return refpTypeName;
898 }
899 
901 {
902  std::string refpTypeName = "";
903  std::stringstream rpName;
904 
905  if(!gel)
906  {
907  std::cout << "Null geoelement parameter on " << __PRETTY_FUNCTION__ << std::endl;
908  DebugStop();
909  }
910 
911  int nsides = gel->NSides();
912  TPZVec<int> sidesToRefine;
913  bool thereIsAnyNeighbourRefined = TPZRefPatternTools::SidesToRefine(gel, sidesToRefine);
914  if(!thereIsAnyNeighbourRefined)
915  {
916  return refpTypeName;
917  }
918 
919  std::string perfix = gel->TypeName();
920  for(int i = 0; i < 3; i++)
921  {
922  rpName << perfix[i];
923  }
924 
925  for(int s = 0; s < nsides; s++)
926  {
927  rpName << sidesToRefine[s];
928  }
929  rpName >> refpTypeName;
930 
931  if(refpTypeName.length() == 0)
932  {
933  refpTypeName = TPZRefPattern::fNonInitializedName;
934  }
935 
936  return refpTypeName;
937 }
938 
940 {
941  bool thereIsAnyNeighbourRefined = false;
942 
943  int ncorners = gel->NCornerNodes();
944  int nsides = gel->NSides();
945  sidestoref.Resize(nsides, 0);
946 
947  for(int s = ncorners; s < nsides; s++)
948  {
949  TPZGeoElSide gelside(gel, s);
950  TPZGeoElSide neighside = gelside.Neighbour();
951  if(!neighside.Exists())
952  {
953  break;
954  }
955  while(neighside != gelside)
956  {
957  if(neighside.Element()->HasSubElement() && neighside.Element()->NSideSubElements(neighside.Side()) > 1)
958  {
959  thereIsAnyNeighbourRefined = true;
960 
961  int ns = neighside.Side();
962  TPZVec<int> MidNodesIndexes;
963 
964  TPZAutoPointer<TPZRefPattern> elrefpattern = neighside.Element()->GetRefPattern();
965  if (!elrefpattern) {
966  DebugStop();
967  }
968  TPZAutoPointer<TPZRefPattern> refSide = elrefpattern->SideRefPattern(ns);
969  if (!refSide)
970  {
971  std::cout << "An element is refined but the pattern has no side refpattern along side " << ns << "\n";
972  neighside.Element()->GetRefPattern()->Print(cout);
973 
974  DebugStop();
975  }
976  elrefpattern->SideNodes(ns, MidNodesIndexes);
977 
978  sidestoref[s] = MidNodesIndexes.NElements();
979 
980  break;
981  }
982  neighside = neighside.Neighbour();
983  }
984  }
985 
986  return thereIsAnyNeighbourRefined;
987 }
988 
989 void TPZRefPatternTools::RefineDirectional(TPZGeoEl *gel, std::set<int> &matids)
990 {
991  if(gel->HasSubElement())
992  {
993  return;
994  }
995  int matid = gel->MaterialId();
996  if(matids.count(matid)) return;
997  TPZManVector<int,27> sidestorefine(gel->NSides(), 0);
998  TPZManVector<int,27> cornerstorefine(gel->NSides(), 0);
999 
1000  //Look for corners which are on the boundary
1001  int numrefribs = 0;
1002  int numrefcorners = 0;
1003  for(int in=0; in<gel->NCornerNodes(); in++)
1004  {
1005  TPZGeoElSide gels(gel,in);
1006  TPZGeoElSide neigh(gels.Neighbour());
1007  while(gels != neigh)
1008  {
1009  if(matids.count(neigh.Element()->MaterialId()))
1010  {
1011  numrefcorners++;
1012  cornerstorefine[in] = 1;
1013  break;
1014  }
1015  neigh = neigh.Neighbour();
1016  }
1017  }
1018 
1019  // nothing to be done
1020  if (numrefcorners == 0) {
1021  return;
1022  }
1023 
1024  // look for ribs which touch the boundary but which do no lay on the boundary
1025  for(int is=gel->NCornerNodes(); is<gel->NSides(); is++)
1026  {
1027  // we are only interested in ribs
1028  if(gel->SideDimension(is) != 1)
1029  {
1030  continue;
1031  }
1032 
1033  // the side is a candidate if it contains a corner which is neighbour of the boundary condition
1034  if(cornerstorefine[gel->SideNodeLocIndex(is,0)] || cornerstorefine[gel->SideNodeLocIndex(is,1)])
1035  {
1036  sidestorefine[is] = 1;
1037  numrefribs++;
1038  TPZGeoElSide gels(gel,is);
1039  TPZGeoElSide neigh(gels.Neighbour());
1040  while(neigh != gels)
1041  {
1042  // if this condition is true the rib lies on the boundary
1043  if(matids.count(neigh.Element()->MaterialId()))
1044  {
1045  sidestorefine[is] = 0;
1046  numrefribs--;
1047  break;
1048  }
1049  neigh = neigh.Neighbour();
1050  }
1051  }
1052  }
1053 
1055  if(!numrefribs)
1056  {
1057  if (gel->Type() == EQuadrilateral && numrefcorners == 4) {
1058  patt = gRefDBase.FindRefPattern(200);
1059 
1060  }
1061  else
1062  {
1063  std::cout << "ncorners to refine " << numrefcorners << " numrefribs " << numrefribs << " i dont understand\n";
1064  return;
1065  }
1066  }
1067  else
1068  {
1069  patt = TPZRefPatternTools::PerfectMatchRefPattern(gel, sidestorefine);
1070  }
1071 
1072  if(patt)
1073  {
1074  gel->SetRefPattern(patt);
1076  gel->Divide(subel);
1077  }
1078  else
1079  {
1080 #ifdef WIN32
1081  DebugStop();
1082 #endif
1083  std::cout << "|"; std::cout.flush();
1084  std::ofstream arquivo ("NotListedPatterns.txt",std::ios::app);
1085  std::list<TPZAutoPointer<TPZRefPattern> >::iterator it;
1086  arquivo << "Compatible refinement patterns\n";
1087 
1088  arquivo << std::endl;
1089  arquivo << "Element Type : " << gel->Type() << ' ' << gel->TypeName() << std::endl;
1090  arquivo << "Sides selected for refinement :" << std::endl;
1091  int i;
1092  for (i=0 ; i<gel->NSides() ; i++)
1093  {
1094  /*
1095  if(cornerstorefine[i] == 1)
1096  {
1097  arquivo << " " << i << " ";
1098  }
1099  */
1100  if (sidestorefine[i] == 1) {
1101  arquivo << " " << i << " " ;
1102  }
1103  }
1104  gel->Print(arquivo);
1105  int in;
1106  arquivo << std::endl;
1107  arquivo << "Neighbouring information \n";
1108  for(in=0; in<gel->NSides(); in++)
1109  {
1110  arquivo << "Side : " << in << " ";
1111  TPZGeoElSide gels(gel,in);
1112  arquivo << "Dim " << gels.Dimension() << " ";
1113  TPZGeoElSide neigh(gels.Neighbour());
1114  while(gels != neigh)
1115  {
1116  if(matids.count(neigh.Element()->MaterialId()))
1117  {
1118  arquivo << neigh.Element()->Id() << "-l-" << neigh.Side() << " ";
1119  if (neigh.Side() == 9 && gel->Type() == ETetraedro)
1120  {
1121  arquivo << "Teje pego meliante..." << std::endl;
1122  neigh.Element()->Print(arquivo);
1123  }
1124  }
1125  neigh = neigh.Neighbour();
1126  }
1127  arquivo << std::endl;
1128  }
1129  arquivo << std::endl;
1130 
1131  arquivo << "Element information : " << gel->Index() << std::endl;
1132  arquivo << "Neighbours of the sides selected for refinement:" << std::endl;
1133  for (i=0 ; i<gel->NSides() ; i++)
1134  {
1135  if(cornerstorefine[i] == 1 || sidestorefine[i] == 1)
1136  {
1137  TPZGeoElSide gelside (gel,i);
1138  TPZGeoElSide neigh = gelside.Neighbour();
1139  while (neigh != gelside)
1140  {
1141  arquivo << "*********** my side = " << i << " neighside " << neigh.Side() << std::endl;
1142  neigh.Element()->Print(arquivo);
1143  neigh = neigh.Neighbour();
1144  }
1145  }
1146  }
1147  arquivo << std::endl << std::endl << std::endl << std::endl;
1148  }
1149 
1150  return;
1151 }
1152 
1153 void TPZRefPatternTools::RefineDirectional(TPZGeoMesh *gmesh, std::set<int> &matids)
1154 {
1155  int64_t nel = gmesh->NElements();
1156  for (int64_t el=0; el<nel; el++) {
1157  TPZGeoEl *gel = gmesh->Element(el);
1158  if (!gel) {
1159  continue;
1160  }
1161  RefineDirectional(gel, matids);
1162  }
1163 }
1164 
1165 void TPZRefPatternTools::RefineDirectional(TPZGeoMesh *gmesh, std::set<int> &matids, int gelMat)
1166 {
1167  int64_t nel = gmesh->NElements();
1168  for (int64_t el=0; el<nel; el++) {
1169  TPZGeoEl *gel = gmesh->Element(el);
1170  if (!gel) {
1171  continue;
1172  }
1173  RefineDirectional(gel, matids, gelMat);
1174  }
1175 }
1176 
1177 
1178 
1179 void TPZRefPatternTools::RefineDirectional(TPZGeoEl *gel, std::set<int> &matids, int gelMat)
1180 {
1181  if(gel->HasSubElement())
1182  {
1183  return;
1184  }
1185  int matid = gel->MaterialId();
1186  if(matids.count(matid)) return;
1187  TPZManVector<int,27> sidestorefine(gel->NSides(), 0);
1188  TPZManVector<int,27> cornerstorefine(gel->NSides(), 0);
1189 
1190  //Look for corners which are on the boundary
1191  int numrefribs = 0;
1192  for(int in=0; in<gel->NCornerNodes(); in++)
1193  {
1194  TPZGeoElSide gels(gel,in);
1195  TPZGeoElSide neigh(gels.Neighbour());
1196  while(gels != neigh)
1197  {
1198  if(matids.count(neigh.Element()->MaterialId()))
1199  {
1200  cornerstorefine[in] = 1;
1201  break;
1202  }
1203  neigh = neigh.Neighbour();
1204  }
1205  }
1206  // look for ribs which touch the boundary but which do not lay on the boundary
1207  for(int is = gel->NCornerNodes(); is < gel->NSides(); is++)
1208  {
1209  // we are only interested in ribs
1210  if(gel->SideDimension(is) != 1)
1211  {
1212  continue;
1213  }
1214 
1215  // the side is a candidate if it contains a corner which is neighbour of the boundary condition
1216  if(cornerstorefine[gel->SideNodeLocIndex(is,0)] || cornerstorefine[gel->SideNodeLocIndex(is,1)])
1217  {
1218  sidestorefine[is] = 1;
1219  numrefribs++;
1220  TPZGeoElSide gels(gel,is);
1221  TPZGeoElSide neigh(gels.Neighbour());
1222  while(neigh != gels)
1223  {
1224  // if this condition is true the rib lies on the boundary
1225  if(matids.count(neigh.Element()->MaterialId()))
1226  {
1227  sidestorefine[is] = 0;
1228  numrefribs--;
1229  break;
1230  }
1231  neigh = neigh.Neighbour();
1232  }
1233  }
1234  }
1235  if(!numrefribs)
1236  {
1237  return;
1238  }
1239 
1241  if(patt)
1242  {
1243  gel->SetMaterialId(gelMat+1);
1244  gel->SetRefPattern(patt);
1246  gel->Divide(subel);
1247  gel->SetMaterialId(gelMat);
1248  }
1249  else
1250  {
1251  std::cout << "|"; std::cout.flush();
1252  static int nTimes {-1};
1253  nTimes++;
1254  std::ofstream arquivo ("NotListedPatterns"+std::to_string(nTimes)+".txt",std::ios::trunc);
1255  std::list<TPZAutoPointer<TPZRefPattern> >::iterator it;
1256  arquivo << "Compatible refinement patterns\n";
1257 
1258  arquivo << std::endl;
1259  arquivo << "Element Type : " << gel->Type() << ' ' << gel->TypeName() << std::endl;
1260  arquivo << "Sides selected for refinement :" << std::endl;
1261  int i;
1262  for (i=0 ; i<gel->NSides() ; i++)
1263  {
1264  if(cornerstorefine[i] == 1)
1265  {
1266  arquivo << " " << i << " ";
1267  }
1268  if (sidestorefine[i] == 1) {
1269  arquivo << " " << i << " " ;
1270  }
1271  }
1272  gel->Print(arquivo);
1273  int in;
1274  arquivo << std::endl;
1275 
1276  arquivo << "Element information : " << gel->Index() << std::endl;
1277  arquivo << "Neighbours of the sides selected for refinement:" << std::endl;
1278  for (i=0 ; i<gel->NSides() ; i++)
1279  {
1280  if(cornerstorefine[i] == 1 || sidestorefine[i] == 1)
1281  {
1282  TPZGeoElSide gelside (gel,i);
1283  TPZGeoElSide neigh = gelside.Neighbour();
1284  while (neigh != gelside)
1285  {
1286  arquivo << "*********** my side = " << i << " neighside " << neigh.Side() << std::endl;
1287  neigh.Element()->Print(arquivo);
1288  neigh = neigh.Neighbour();
1289  }
1290  }
1291  }
1292  arquivo << std::endl << std::endl << std::endl << std::endl;
1293  }
1294 
1295  return;
1296 }
1297 
1299 {
1300  int nsides = gel->NSides();
1301  for(int s = 0; s < nsides; s++)
1302  {
1303  TPZGeoElSide gelside(gel,s);
1304  gelside = gelside.LowestFatherSide();
1305  TPZGeoElSide neighside(gelside.Neighbour());
1306  while(gelside != neighside)
1307  {
1308  if(matids.count(neighside.Element()->MaterialId()))
1309  {
1311  if(refP)
1312  {
1313  TPZVec<TPZGeoEl*> sons;
1314  gel->SetRefPattern(refP);
1315  gel->Divide(sons);
1316 
1317  return;
1318  }
1319  else
1320  {
1321  std::cout << "Uniform refpattern was not found!" << std::endl;
1322  std::cout << "See " << __PRETTY_FUNCTION__ << std::endl;
1323  return;
1324  }
1325  }
1326  neighside = neighside.Neighbour();
1327  }
1328  }
1329 }
1330 
1332 {
1333  TPZIntPoints * rule = gelside.Element()->CreateSideIntegrationRule(gelside.Side(),5);
1334  int dim = gelside.Dimension();
1335  TPZManVector<REAL,4> pt(dim,0.);
1336  REAL weight;
1337  int np = rule->NPoints();
1338  rule->Point(0, pt, weight);
1339  TPZFNMatrix<9> jacobian(dim,dim,0.),jacinv(dim,dim,0.);
1340  TPZFNMatrix<9> axes(3,3,0.);
1341  REAL detjac;
1342  gelside.Jacobian(pt, jacobian, axes , detjac, jacinv);
1343  int ip;
1344  for(ip=1; ip<np; ip++)
1345  {
1346  rule->Point(ip, pt, weight);
1347  TPZFNMatrix<9> jacobian2(dim,dim,0.),jacinv2(dim,dim,0.);
1348  REAL detjac2;
1349  TPZFNMatrix<9> axes2(3,3,0.);
1350  gelside.Jacobian(pt, jacobian2, axes2 , detjac2, jacinv2);
1351  jacobian2 -= jacobian;
1352  axes2 -= axes;
1353  detjac2 -= detjac;
1354  jacinv2 -= jacinv;
1355  REAL diff = Norm(jacobian2)+Norm(axes2)+Norm(jacinv2)+fabs(detjac2);
1356  if(diff > 1.e-8)
1357  {
1358  delete rule;
1359  return false;
1360  }
1361  }
1362  delete rule;
1363 
1364  return true;
1365 }
1366 
1368 {
1369  int isub,sd,ip;
1370 
1371  //x1 no filho deformado, x2 no pai deformado
1372  TPZManVector<REAL,3> x1(3),pf(3),x2(3),xpf(3,0.);
1373  REAL weight;
1374  TPZGeoEl *father = refp->Element(0);//pai
1375  int dimfatside,fatside,nsides;
1376  int dimfat = father->Dimension();
1377  REAL Tol;
1378  ZeroTolerance(Tol);
1379  if(Tol < 1.e-10) Tol = 1.e-10;
1380 
1381  TPZGeoEl *subel;
1382  int nsubs = refp->NSubElements();
1383  for(isub=0;isub<nsubs;isub++)
1384  {
1385  subel = refp->Element(isub+1);
1386  nsides = subel->NSides();
1387  for(sd=0; sd<nsides; sd++)
1388  {
1389  TPZGeoElSide subside(subel,sd);
1390  if(!ConstJacobian(subside,1.e-6))
1391  {
1392  continue;
1393  }
1394  if( sd<subel->NNodes() && refp->IsFatherNeighbour(TPZGeoElSide(subel,sd),father) ) continue;
1395 
1396  int dims = subel->SideDimension(sd);
1397  int dimsub = subel->Dimension();
1398 
1399  //regra de integracao para o espaco paramcorico do lado do sub-elemento
1400  TPZIntPoints *rule = subel->CreateSideIntegrationRule(sd,5);
1401  TPZVec<int> order(dims,2);
1402  TPZManVector<REAL,3> point(dims,0.),point2(dimsub),pointparamfather(dimfat,0.);
1403  rule->SetOrder(order);
1404 
1405  for(ip=0;ip<rule->NPoints();ip++)
1406  {
1407  //ponto no espaco paramcorico do lado do filho
1408  rule->Point(ip,point,weight);
1409  TPZTransform<> sidet(dimsub);//transformacao unitcoia
1410  if(dims < dimsub)
1411  {
1412  sidet = subel->SideToSideTransform(sd,nsides-1);//transf. no elemento metre
1413  }
1414 
1415  //transformacao para o interior do mestre do sub-elemento
1416  sidet.Apply(point,point2);
1417 
1418  //ponto no lado do filho deformado
1419  subel->X(point2,x1);
1420 
1421  //transformacao: espaco paramcorico do filho/lado -> espaco paramcorico do pai/lado
1422  father->ComputeXInverse(x1, pointparamfather, Tol);
1423 
1424  //no elemento mestre do pai, father coo deformado
1425  fatside = father->WhichSide(pointparamfather);
1426  dimfatside = father->SideDimension(fatside);
1427 
1428  //transformacoo calculada pelo TPZRefPattern
1429  TPZTransform<> trans = refp->Transform(sd,isub);
1430 
1431  TPZVec<REAL> pf(dimfatside);
1432 
1433  //ponto pf no espaco parametrico do lado do pai
1434  trans.Apply(point,pf);
1435 
1436  //unitcoia do lado do pai
1437  TPZTransform<> sidetf(dimfatside);
1438  if(dimfatside < dimfat)
1439  {
1440  //para o interior do sub-elemento
1441  sidetf = father->SideToSideTransform(fatside,father->NSides()-1);
1442 
1443  //sidetf = father->SideToElemShapeT(fatside);
1444  }
1445 
1446  //do espaco parametrico do lado do pai para o interior do pai
1447  sidetf.Apply(pf,xpf);
1448 
1449  //ponto no lado do pai deformado
1450  father->X(xpf,x2);
1451  if( sqrt( (x1[0]-x2[0])*(x1[0]-x2[0]) + (x1[1]-x2[1])*(x1[1]-x2[1]) + (x1[2]-x2[2])*(x1[2]-x2[2]) ) > 1.e-10 )
1452  {
1453  PZError << "\nTransformacao errada\n";
1454  PZError << "son = " << (subel->Id()) << std::endl;
1455  PZError << "father = " << (father->Id()) << std::endl;
1456  PZError << "side = " << sd << std::endl << std::endl;
1457 
1458  DebugStop();
1459  }
1460  else
1461  {
1462  // std::cout << "Transformacao OK!\n";
1463  // std::cout << "Filho/lado : " << subel->Id() << "/" << sd << std::endl;
1464  // std::cout << "Pai : " << father->Id() << std::endl << std::endl;
1465  }//fim if sqrt..
1466  }//fim rule
1467 
1468  delete rule;
1469  }//fim for sd
1470  }//fim for isub
1471 }
1472 
1473 void TPZRefPatternTools::NodesHunter(TPZGeoMesh &gMesh, TPZVec<int>& NodesHunted, int IdIni, int IdFin, double Tol)
1474 {
1477  int dim = 3;
1478 
1479  int VecSize = gMesh.NodeVec().NElements();
1480  int posIni = -1;
1481  int posFin = -1;
1482 
1485  TPZFMatrix<REAL> VectorialNotation(VecSize,dim);
1486  for(int pos = 0; pos < VecSize; pos++)
1487  {
1488  if(gMesh.NodeVec()[pos].Id() == IdIni) posIni = pos;
1489  if(gMesh.NodeVec()[pos].Id() == IdFin) posFin = pos;
1490 
1491  for(int coord = 0; coord < dim; coord++)
1492  {
1493  double val = gMesh.NodeVec()[pos].Coord(coord) - gMesh.NodeVec()[IdIni].Coord(coord);
1494  VectorialNotation.PutVal(pos,coord,val);
1495  }
1496  }
1497 
1498 #ifdef PZDEBUG
1499  if(posIni == -1 || posFin == -1)
1500  {
1501  std::cout << "Initial Node index or Final Node index doesn't belong to the given TPZGeoNode TPZVec!\n";
1502  std::cout << "See NodesHunter method!\n";
1503  DebugStop();
1504  }
1505 #endif
1506 
1507  // Computing BasisChange Matrix
1508  // Where NewBase X_axis is defined by IniNode->FinNode orientation
1509  // and NewBase Y_axis is perpendicular to X_axis in XY plane counter-clockwise
1510  TPZFMatrix<REAL> IfromCntoBase(dim,dim,0.);
1511  IfromCntoBase(dim-1, dim-1) = 1.;
1512  double norm = 0.;
1513 
1514  //computing X_axis
1515  for(int coord = 0; coord < dim; coord++)
1516  {
1517  double val = VectorialNotation.GetVal(posFin,coord) - VectorialNotation.GetVal(posIni,coord);
1518  IfromCntoBase.PutVal(0,coord,val);
1519  norm += val*val;
1520  }
1521  //normalizing X_axis
1522  for(int coord = 0; coord < dim; coord++)
1523  {
1524  IfromCntoBase.PutVal(0,coord,IfromCntoBase.GetVal(0,coord)/sqrt(norm));
1525  }
1526  //computing Y_axis, i.e., if X_axis=(x,y) then Y_axis=(-y,x)
1527  IfromCntoBase.PutVal(1,0,-IfromCntoBase.GetVal(0,1));
1528  IfromCntoBase.PutVal(1,1, IfromCntoBase.GetVal(0,0));
1529 
1531  for(int i = 0; i < VecSize; i++)
1532  {
1533  TPZManVector <double> temp; temp.Resize(3);
1534  for(int j = 0; j < dim; j++)
1535  {
1536  double val = 0.;
1537  for(int k = 0; k < dim; k++)
1538  {
1539  val += IfromCntoBase.GetVal(j,k)*VectorialNotation.GetVal(i,k);
1540  }
1541  temp[j] = val;
1542  }
1543  for(int p = 0; p < dim; p++) VectorialNotation.PutVal(i,p,temp[p]);
1544  }
1545 
1547  std::map <double,int> mymap;
1548  for(int h = 0; h < VecSize; h++)
1549  {
1550  if( VectorialNotation.GetVal(h,0) >= VectorialNotation.GetVal(posIni,0) &&
1551  VectorialNotation.GetVal(h,0) <= VectorialNotation.GetVal(posFin,0) &&
1552  (fabs(VectorialNotation.GetVal(h,1)) + fabs(VectorialNotation.GetVal(h,2))) < Tol )
1553  {
1554  std::pair< double , int> Item(VectorialNotation.GetVal(h,0), gMesh.NodeVec()[h].Id());
1555  mymap.insert(Item);
1556  }
1557  }
1558  NodesHunted.Resize(mymap.size());
1559  std::map<double, int>::iterator it = mymap.begin();
1560  int i = 0;
1561  for(it = mymap.begin(); it!= mymap.end(); it++, i++)
1562  {
1563  NodesHunted[i] = it->second;
1564  }
1565 }
1566 
1568 {
1569  int id;
1570  GetElTypePermutations(gel->Type(), permutation);
1571  for(int p = 0; p < permutation.NElements(); p++)
1572  {
1573  for(int n = 0; n < permutation[p].NElements(); n++)
1574  {
1575  id = gel->NodeIndex(permutation[p][n]);
1576  permutation[p][n] = id;
1577  }
1578  }
1579 }
1580 
1582 {
1583  int nperm, nnodes, p;
1584  switch (elType)
1585  {
1586  case (EOned) :
1587  {
1588  nperm = 2;
1589  nnodes = 2;
1590  permutation.Resize(nperm);
1591  for(p = 0; p < nperm; p++)
1592  {
1593  permutation[p].Resize(nnodes);
1594  }
1595  p = 0;
1596  permutation[p][0] = 0; permutation[p][1] = 1; p++;
1597  permutation[p][0] = 1; permutation[p][1] = 0;
1598 
1599  break;
1600  }
1601  case (ETriangle) :
1602  {
1603  nperm = 6;
1604  nnodes = 3;
1605  permutation.Resize(nperm);
1606  for(p = 0; p < nperm; p++)
1607  {
1608  permutation[p].Resize(nnodes);
1609  }
1610  p = 0;
1611  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; p++;
1612  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 0; p++;
1613  permutation[p][0] = 2; permutation[p][1] = 0; permutation[p][2] = 1; p++;
1614 
1615  permutation[p][0] = 0; permutation[p][1] = 2; permutation[p][2] = 1; p++;
1616  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; p++;
1617  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 2;
1618 
1619  break;
1620  }
1621  case (EQuadrilateral) :
1622  {
1623  nperm = 8;
1624  nnodes = 4;
1625  permutation.Resize(nperm);
1626  for(p = 0; p < nperm; p++)
1627  {
1628  permutation[p].Resize(nnodes);
1629  }
1630  p = 0;
1631  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 3; p++;
1632  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 0; p++;
1633  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 1; p++;
1634  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 2; p++;
1635 
1636  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 1; p++;
1637  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 0; p++;
1638  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 3; p++;
1639  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 2;
1640 
1641  break;
1642  }
1643  case (ETetraedro) :
1644  {
1645  nperm = 24;
1646  nnodes = 4;
1647  permutation.Resize(nperm);
1648  for(p = 0; p < nperm; p++)
1649  {
1650  permutation[p].Resize(nnodes);
1651  }
1652  p = 0;
1653  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 3; p++;
1654  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 0; permutation[p][3] = 3; p++;
1655  permutation[p][0] = 2; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 3; p++;
1656 
1657  permutation[p][0] = 0; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 3; p++;
1658  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 3; p++;
1659  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 2; permutation[p][3] = 3; p++;
1660  //
1661  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 3; permutation[p][3] = 2; p++;
1662  permutation[p][0] = 1; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 2; p++;
1663  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 2; p++;
1664 
1665  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 1; permutation[p][3] = 2; p++;
1666  permutation[p][0] = 3; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 2; p++;
1667  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 2; p++;
1668  //
1669  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 1; p++;
1670  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 0; permutation[p][3] = 1; p++;
1671  permutation[p][0] = 2; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 1; p++;
1672 
1673  permutation[p][0] = 0; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 1; p++;
1674  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 1; p++;
1675  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 2; permutation[p][3] = 1; p++;
1676  //
1677  permutation[p][0] = 3; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 0; p++;
1678  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 0; p++;
1679  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 1; permutation[p][3] = 0; p++;
1680 
1681  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 0; p++;
1682  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 3; permutation[p][3] = 0; p++;
1683  permutation[p][0] = 1; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 0;
1684 
1685  break;
1686  }
1687  case (EPiramide) :
1688  {
1689  nperm = 8;
1690  nnodes = 5;
1691  permutation.Resize(nperm);
1692  for(p = 0; p < nperm; p++)
1693  {
1694  permutation[p].Resize(nnodes);
1695  }
1696  p = 0;
1697  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 3; permutation[p][4] = 4; p++;
1698  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 0; permutation[p][4] = 4; p++;
1699  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 1; permutation[p][4] = 4; p++;
1700  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 2; permutation[p][4] = 4; p++;
1701 
1702  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 1; permutation[p][4] = 4; p++;
1703  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 0; permutation[p][4] = 4; p++;
1704  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 3; permutation[p][4] = 4; p++;
1705  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 2; permutation[p][4] = 4;
1706 
1707  break;
1708  }
1709  case (EPrisma) :
1710  {
1711  nperm = 12;
1712  nnodes = 6;
1713  permutation.Resize(nperm);
1714  for(p = 0; p < nperm; p++)
1715  {
1716  permutation[p].Resize(nnodes);
1717  }
1718  p = 0;
1719  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 3; permutation[p][4] = 4; permutation[p][5] = 5; p++;
1720  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 0; permutation[p][3] = 4; permutation[p][4] = 5; permutation[p][5] = 3; p++;
1721  permutation[p][0] = 2; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 5; permutation[p][4] = 3; permutation[p][5] = 4; p++;
1722 
1723  permutation[p][0] = 0; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 3; permutation[p][4] = 5; permutation[p][5] = 4; p++;
1724  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 5; permutation[p][4] = 4; permutation[p][5] = 3; p++;
1725  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 2; permutation[p][3] = 4; permutation[p][4] = 3; permutation[p][5] = 5; p++;
1726  //
1727  permutation[p][0] = 3; permutation[p][1] = 4; permutation[p][2] = 5; permutation[p][3] = 0; permutation[p][4] = 1; permutation[p][5] = 2; p++;
1728  permutation[p][0] = 4; permutation[p][1] = 5; permutation[p][2] = 3; permutation[p][3] = 1; permutation[p][4] = 2; permutation[p][5] = 0; p++;
1729  permutation[p][0] = 5; permutation[p][1] = 3; permutation[p][2] = 4; permutation[p][3] = 2; permutation[p][4] = 0; permutation[p][5] = 1; p++;
1730 
1731  permutation[p][0] = 3; permutation[p][1] = 5; permutation[p][2] = 4; permutation[p][3] = 0; permutation[p][4] = 2; permutation[p][5] = 1; p++;
1732  permutation[p][0] = 5; permutation[p][1] = 4; permutation[p][2] = 3; permutation[p][3] = 2; permutation[p][4] = 1; permutation[p][5] = 0; p++;
1733  permutation[p][0] = 4; permutation[p][1] = 3; permutation[p][2] = 5; permutation[p][3] = 1; permutation[p][4] = 0; permutation[p][5] = 2;
1734 
1735  break;
1736  }
1737  case (ECube) :
1738  {
1739  nperm = 48;
1740  nnodes = 8;
1741  permutation.Resize(nperm);
1742  for(p = 0; p < nperm; p++)
1743  {
1744  permutation[p].Resize(nnodes);
1745  }
1746  p = 0;
1747  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 3; permutation[p][4] = 4; permutation[p][5] = 5; permutation[p][6] = 6; permutation[p][7] = 7; p++;
1748  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 0; permutation[p][4] = 5; permutation[p][5] = 6; permutation[p][6] = 7; permutation[p][7] = 4; p++;
1749  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 1; permutation[p][4] = 6; permutation[p][5] = 7; permutation[p][6] = 4; permutation[p][7] = 5; p++;
1750  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 2; permutation[p][4] = 7; permutation[p][5] = 4; permutation[p][6] = 5; permutation[p][7] = 6; p++;
1751 
1752  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 1; permutation[p][4] = 4; permutation[p][5] = 7; permutation[p][6] = 6; permutation[p][7] = 5; p++;
1753  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 0; permutation[p][4] = 7; permutation[p][5] = 6; permutation[p][6] = 5; permutation[p][7] = 4; p++;
1754  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 3; permutation[p][4] = 6; permutation[p][5] = 5; permutation[p][6] = 4; permutation[p][7] = 7; p++;
1755  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 2; permutation[p][4] = 5; permutation[p][5] = 4; permutation[p][6] = 7; permutation[p][7] = 6; p++;
1756  //
1757  permutation[p][0] = 0; permutation[p][1] = 1; permutation[p][2] = 5; permutation[p][3] = 4; permutation[p][4] = 3; permutation[p][5] = 2; permutation[p][6] = 6; permutation[p][7] = 7; p++;
1758  permutation[p][0] = 1; permutation[p][1] = 5; permutation[p][2] = 4; permutation[p][3] = 0; permutation[p][4] = 2; permutation[p][5] = 6; permutation[p][6] = 7; permutation[p][7] = 3; p++;
1759  permutation[p][0] = 5; permutation[p][1] = 4; permutation[p][2] = 0; permutation[p][3] = 1; permutation[p][4] = 6; permutation[p][5] = 7; permutation[p][6] = 3; permutation[p][7] = 2; p++;
1760  permutation[p][0] = 4; permutation[p][1] = 0; permutation[p][2] = 1; permutation[p][3] = 5; permutation[p][4] = 7; permutation[p][5] = 3; permutation[p][6] = 2; permutation[p][7] = 6; p++;
1761 
1762  permutation[p][0] = 0; permutation[p][1] = 4; permutation[p][2] = 5; permutation[p][3] = 1; permutation[p][4] = 3; permutation[p][5] = 7; permutation[p][6] = 6; permutation[p][7] = 2; p++;
1763  permutation[p][0] = 4; permutation[p][1] = 5; permutation[p][2] = 1; permutation[p][3] = 0; permutation[p][4] = 7; permutation[p][5] = 6; permutation[p][6] = 2; permutation[p][7] = 3; p++;
1764  permutation[p][0] = 5; permutation[p][1] = 1; permutation[p][2] = 0; permutation[p][3] = 4; permutation[p][4] = 6; permutation[p][5] = 2; permutation[p][6] = 3; permutation[p][7] = 7; p++;
1765  permutation[p][0] = 1; permutation[p][1] = 0; permutation[p][2] = 4; permutation[p][3] = 5; permutation[p][4] = 2; permutation[p][5] = 3; permutation[p][6] = 7; permutation[p][7] = 6; p++;
1766  //
1767  permutation[p][0] = 1; permutation[p][1] = 2; permutation[p][2] = 6; permutation[p][3] = 5; permutation[p][4] = 0; permutation[p][5] = 3; permutation[p][6] = 7; permutation[p][7] = 4; p++;
1768  permutation[p][0] = 2; permutation[p][1] = 6; permutation[p][2] = 5; permutation[p][3] = 1; permutation[p][4] = 3; permutation[p][5] = 7; permutation[p][6] = 4; permutation[p][7] = 0; p++;
1769  permutation[p][0] = 6; permutation[p][1] = 5; permutation[p][2] = 1; permutation[p][3] = 2; permutation[p][4] = 7; permutation[p][5] = 4; permutation[p][6] = 0; permutation[p][7] = 3; p++;
1770  permutation[p][0] = 5; permutation[p][1] = 1; permutation[p][2] = 2; permutation[p][3] = 6; permutation[p][4] = 4; permutation[p][5] = 0; permutation[p][6] = 3; permutation[p][7] = 7; p++;
1771 
1772  permutation[p][0] = 1; permutation[p][1] = 5; permutation[p][2] = 6; permutation[p][3] = 2; permutation[p][4] = 0; permutation[p][5] = 4; permutation[p][6] = 7; permutation[p][7] = 3; p++;
1773  permutation[p][0] = 5; permutation[p][1] = 6; permutation[p][2] = 2; permutation[p][3] = 1; permutation[p][4] = 4; permutation[p][5] = 7; permutation[p][6] = 3; permutation[p][7] = 0; p++;
1774  permutation[p][0] = 6; permutation[p][1] = 2; permutation[p][2] = 1; permutation[p][3] = 5; permutation[p][4] = 7; permutation[p][5] = 3; permutation[p][6] = 0; permutation[p][7] = 4; p++;
1775  permutation[p][0] = 2; permutation[p][1] = 1; permutation[p][2] = 5; permutation[p][3] = 6; permutation[p][4] = 3; permutation[p][5] = 0; permutation[p][6] = 4; permutation[p][7] = 7; p++;
1776  //
1777  permutation[p][0] = 2; permutation[p][1] = 3; permutation[p][2] = 7; permutation[p][3] = 6; permutation[p][4] = 1; permutation[p][5] = 0; permutation[p][6] = 4; permutation[p][7] = 5; p++;
1778  permutation[p][0] = 3; permutation[p][1] = 7; permutation[p][2] = 6; permutation[p][3] = 2; permutation[p][4] = 0; permutation[p][5] = 4; permutation[p][6] = 5; permutation[p][7] = 1; p++;
1779  permutation[p][0] = 7; permutation[p][1] = 6; permutation[p][2] = 2; permutation[p][3] = 3; permutation[p][4] = 4; permutation[p][5] = 5; permutation[p][6] = 1; permutation[p][7] = 0; p++;
1780  permutation[p][0] = 6; permutation[p][1] = 2; permutation[p][2] = 3; permutation[p][3] = 7; permutation[p][4] = 5; permutation[p][5] = 1; permutation[p][6] = 0; permutation[p][7] = 4; p++;
1781 
1782  permutation[p][0] = 2; permutation[p][1] = 6; permutation[p][2] = 7; permutation[p][3] = 3; permutation[p][4] = 1; permutation[p][5] = 5; permutation[p][6] = 4; permutation[p][7] = 0; p++;
1783  permutation[p][0] = 6; permutation[p][1] = 7; permutation[p][2] = 3; permutation[p][3] = 2; permutation[p][4] = 5; permutation[p][5] = 4; permutation[p][6] = 0; permutation[p][7] = 1; p++;
1784  permutation[p][0] = 7; permutation[p][1] = 3; permutation[p][2] = 2; permutation[p][3] = 6; permutation[p][4] = 4; permutation[p][5] = 0; permutation[p][6] = 1; permutation[p][7] = 5; p++;
1785  permutation[p][0] = 3; permutation[p][1] = 2; permutation[p][2] = 6; permutation[p][3] = 7; permutation[p][4] = 0; permutation[p][5] = 1; permutation[p][6] = 5; permutation[p][7] = 4; p++;
1786  //
1787  permutation[p][0] = 3; permutation[p][1] = 0; permutation[p][2] = 4; permutation[p][3] = 7; permutation[p][4] = 1; permutation[p][5] = 5; permutation[p][6] = 6; permutation[p][7] = 4; p++;
1788  permutation[p][0] = 0; permutation[p][1] = 4; permutation[p][2] = 7; permutation[p][3] = 3; permutation[p][4] = 5; permutation[p][5] = 6; permutation[p][6] = 4; permutation[p][7] = 1; p++;
1789  permutation[p][0] = 4; permutation[p][1] = 7; permutation[p][2] = 3; permutation[p][3] = 0; permutation[p][4] = 6; permutation[p][5] = 4; permutation[p][6] = 1; permutation[p][7] = 5; p++;
1790  permutation[p][0] = 7; permutation[p][1] = 3; permutation[p][2] = 0; permutation[p][3] = 4; permutation[p][4] = 4; permutation[p][5] = 1; permutation[p][6] = 5; permutation[p][7] = 6; p++;
1791 
1792  permutation[p][0] = 3; permutation[p][1] = 7; permutation[p][2] = 4; permutation[p][3] = 0; permutation[p][4] = 2; permutation[p][5] = 6; permutation[p][6] = 5; permutation[p][7] = 1; p++;
1793  permutation[p][0] = 7; permutation[p][1] = 4; permutation[p][2] = 0; permutation[p][3] = 3; permutation[p][4] = 6; permutation[p][5] = 5; permutation[p][6] = 1; permutation[p][7] = 2; p++;
1794  permutation[p][0] = 4; permutation[p][1] = 0; permutation[p][2] = 3; permutation[p][3] = 7; permutation[p][4] = 5; permutation[p][5] = 1; permutation[p][6] = 2; permutation[p][7] = 6; p++;
1795  permutation[p][0] = 0; permutation[p][1] = 3; permutation[p][2] = 7; permutation[p][3] = 4; permutation[p][4] = 1; permutation[p][5] = 2; permutation[p][6] = 6; permutation[p][7] = 5; p++;
1796  //
1797  permutation[p][0] = 4; permutation[p][1] = 5; permutation[p][2] = 6; permutation[p][3] = 7; permutation[p][4] = 0; permutation[p][5] = 1; permutation[p][6] = 2; permutation[p][7] = 3; p++;
1798  permutation[p][0] = 5; permutation[p][1] = 6; permutation[p][2] = 7; permutation[p][3] = 4; permutation[p][4] = 1; permutation[p][5] = 2; permutation[p][6] = 3; permutation[p][7] = 0; p++;
1799  permutation[p][0] = 6; permutation[p][1] = 7; permutation[p][2] = 4; permutation[p][3] = 5; permutation[p][4] = 2; permutation[p][5] = 3; permutation[p][6] = 0; permutation[p][7] = 1; p++;
1800  permutation[p][0] = 7; permutation[p][1] = 4; permutation[p][2] = 5; permutation[p][3] = 6; permutation[p][4] = 3; permutation[p][5] = 0; permutation[p][6] = 1; permutation[p][7] = 2; p++;
1801 
1802  permutation[p][0] = 4; permutation[p][1] = 7; permutation[p][2] = 6; permutation[p][3] = 5; permutation[p][4] = 0; permutation[p][5] = 3; permutation[p][6] = 2; permutation[p][7] = 1; p++;
1803  permutation[p][0] = 7; permutation[p][1] = 6; permutation[p][2] = 5; permutation[p][3] = 4; permutation[p][4] = 3; permutation[p][5] = 2; permutation[p][6] = 1; permutation[p][7] = 0; p++;
1804  permutation[p][0] = 6; permutation[p][1] = 5; permutation[p][2] = 4; permutation[p][3] = 7; permutation[p][4] = 2; permutation[p][5] = 1; permutation[p][6] = 0; permutation[p][7] = 3; p++;
1805  permutation[p][0] = 5; permutation[p][1] = 4; permutation[p][2] = 7; permutation[p][3] = 6; permutation[p][4] = 1; permutation[p][5] = 0; permutation[p][6] = 3; permutation[p][7] = 2; p++;
1806 
1807  break;
1808  }
1809  default:
1810  {
1811  cout << "Cant return permutation because MElementType was not found on " << __PRETTY_FUNCTION__ << endl;
1812  DebugStop();
1813  }
1814  }
1815 }
1816 
1817 
1818 
static void NodesHunter(TPZGeoMesh &gMesh, TPZVec< int > &NodesHunted, int IdIni, int IdFin, double Tol=1.E-1)
NodesHunted vector is the sequential nodesIds that belongs (i.e.: "Tol" far) to InitialNode(IdIni)~Fi...
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
virtual int NPoints() const =0
Returns number of points for the cubature rule related.
TPZGeoNode * NodePtr(int i) const
Returns a pointer to the ith node of the element.
Definition: pzgeoel.cpp:2566
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
int MaterialId() const
Returns the material index of the element.
Definition: pzgeoel.h:250
const std::list< TPZAutoPointer< TPZRefPattern > > & RefPatternList(MElementType eltype)
Return the complete set of refinement patterns availabe.
bool IsFatherNeighbour(TPZGeoElSide fatherSide, TPZGeoEl *son) const
Verifies the neighbouring relationship between a son and a father&#39;s side.
int NSubElements()
return the number of element/side pairs which compose the current set of points
MElementType Type()
static std::string BuildRefPatternModelName(TPZRefPattern &refp)
Returns the the name of refpattern model.
void InsertRefPattern(TPZAutoPointer< TPZRefPattern > &refpat)
Insert the refinement pattern in the list of availabe refinement patterns assigns an Id to refPattern...
TPZGeoEl * Element(int iel)
It returns the element number iel from the stack of elements of the geometric mesh.
virtual int NCornerNodes() const =0
Returns the number of corner nodes of the element.
int64_t NElements() const
Number of elements of the mesh.
Definition: pzgmesh.h:129
void InsertPermuted()
Generate all permuted partitions and insert them in the mesh.
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
static void GenerateGMeshFromElementVec(const TPZVec< TPZGeoEl *> &elementVec, TPZGeoMesh &refGMesh)
clarg::argBool h("-h", "help message", false)
int64_t NElements() const
Access method to query the number of elements of the vector.
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZGeoElSide Neighbour() const
Definition: pzgeoel.h:754
TPZGeoMesh * Mesh() const
Returns the mesh to which the element belongs.
Definition: pzgeoel.h:220
virtual int NSides() const =0
Returns the number of connectivities of the element.
TPZTransform Transform(int subElSide, int sub)
It returns the TPZTransform associated with a certain sub-element&#39;s side to the father&#39;s coordinates...
virtual void Divide(TPZVec< TPZGeoEl *> &pv)
Divides the element and puts the resulting elements in the vector.
Definition: pzgeoel.h:742
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...
virtual int SideDimension(int side) const =0
Return the dimension of side.
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZAutoPointer< TPZRefPattern > FindRefPattern(TPZAutoPointer< TPZRefPattern > &refpat)
Check whether the refinement pattern already exists.
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
TPZRefPatternDataBase gRefDBase
External variable to data base of patterns.
static void GetCompatibleRefPatterns(TPZGeoEl *gel, std::list< TPZAutoPointer< TPZRefPattern > > &refs)
Search for refpatterns that could be used by a given element with respect to their neighbours...
TPZGeoNode * FindNode(TPZVec< REAL > &co)
Returns the nearest node to the coordinate. This method is VERY INEFFICIENT.
Definition: pzgmesh.cpp:419
TPZAutoPointer< TPZRefPattern > SideRefPattern(int side)
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
virtual std::string TypeName() const
Returns the type of the element as a string.
Definition: pzgeoel.h:262
virtual TPZIntPoints * CreateSideIntegrationRule(int side, int order)=0
Creates an integration rule for the topology of the corresponding side and able to integrate a polyno...
virtual int NSideSubElements(int side) const =0
Returns the number of subelements as returned by GetSubElements2(side)
TPZGeoElSide LowestFatherSide()
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
static TPZAutoPointer< TPZRefPattern > GetRefPatternBasedOnRealMeshElements(TPZVec< TPZGeoEl *> &realMeshElementVec)
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 Resize(const int newsize)
Increase the size of the chunk vector.
Definition: pzadmchunk.h:280
static void PairMeshesNodesMatchingCoordinates(TPZGeoMesh &meshA, TPZGeoMesh &meshB, TPZTransform<> &fromAtoB, std::map< int, int > &pairedNodes)
This method pair nodes from refA->mesh to refB->mesh, using the givem transformation from refA->mesh ...
static const double tol
Definition: pzgeoprism.cpp:23
TPZAutoPointer< TPZRefPattern > GetUniformRefPattern(MElementType type)
Retrieves the uniform refinement pattern for given element type.
static TPZAutoPointer< TPZRefPattern > DragModelPatNodes(TPZGeoEl *gel, TPZAutoPointer< TPZRefPattern > modelPat, std::map< int, std::pair< TPZGeoEl *, std::map< int, int > > > &neighCorresp)
Return an refpattern based on a gived one (modelPat), whose midnodes was dragged to match with a geoe...
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
TPZGeoEl * Element(int64_t iel)
Definition: pzgmesh.h:139
int64_t Index() const
Returns the index of the element within the element vector of the mesh.
Definition: pzgeoel.h:730
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
int64_t NNodes() const
Number of nodes of the mesh.
Definition: pzgmesh.h:126
static bool SidesToRefine(TPZGeoEl *gel, TPZVec< int > &sidestoref)
Returns if there is any neigbour already refined.
virtual void SetOrder(TPZVec< int > &ord, int type=0)=0
Sets the order of the cubature rule.
TPZTransform< REAL > NeighbourSideTransform(const TPZGeoElSide &neighbour)
virtual MElementType Type() const =0
Returns the element type acording to pzeltype.h.
virtual int SideNodeLocIndex(int side, int nodenum) const =0
Returns the local index of a node on a side.
TPZAdmChunkVector< TPZGeoNode > & NodeVec()
Definition: pzgmesh.h:140
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
void SideNodes(int fatherSide, TPZVec< int > &vecNodes)
It returns a TPZVec containing the nodes contained in a given father&#39;s side.
virtual TPZGeoEl * ClonePatchEl(TPZGeoMesh &DestMesh, std::map< int64_t, int64_t > &gl2lcNdIdx, std::map< int64_t, int64_t > &gl2lcElIdx) const =0
Creates a clone of this element into a new patch mesh.
static void RefineUniformIfNeighMat(TPZGeoEl *gel, std::set< int > &matids)
static void GetElTypePermutations(MElementType elType, TPZVec< TPZManVector< int, 8 > > &permutation)
Fill the TPZVec "permutation" with the valid permutations of a given element type.
Contains the TPZRefPattern class which defines the topology of the current refinement pattern to a me...
int NSideNodes(int fatherSide) const
Returns the number of internal nodes of side.
static bool ConstJacobian(TPZGeoElSide gelside, REAL tol=1.e-6)
Method to test if the jacobian of a TPZGeoElSide element is constant.
virtual int HasSubElement() const =0
Return 1 if the element has subelements.
static void RefineDirectional(TPZGeoEl *gel, std::set< int > &matids)
Refines the element if it touches an element with a material id included in matids.
virtual int NNodes() const =0
Returns the number of nodes of the element.
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
static TPZAutoPointer< TPZRefPattern > PerfectMatchRefPattern(TPZGeoEl *gel, TPZVec< int > &sidestorefine)
Returns the refpattern that matches the sides refinement intensity and midnodes coordinates with resp...
static void PairMeshesCornerNodesMatchingCoordinates(TPZGeoMesh &meshA, TPZGeoMesh &meshB, TPZTransform<> &fromAtoB, std::map< int, int > &pairedNodes)
This method pair CORNER nodes from refA->mesh.father to refB->mesh.father, using the givem transforma...
Contains the TPZRefPatternDataBase class which defines data base of patterns.
virtual int Dimension() const =0
Returns the dimension of the element.
virtual TPZTransform< REAL > SideToSideTransform(int sidefrom, int sideto)=0
Compute the transformation between the master element space of one side of an element to the master e...
Defines the topology of the current refinement pattern to a mesh. Refine.
Definition: TPZRefPattern.h:77
static void GetGelPermutations(TPZGeoEl *gel, TPZVec< TPZManVector< int, 8 > > &permutation)
Fill the TPZVec "permutation" with the valid permutations of "gel".
int64_t Id() const
Returns the Id of the element.
Definition: pzgeoel.h:223
static void PrintGMeshVTKneighbourhood(TPZGeoMesh *gmesh, int64_t elIndex, std::ofstream &file)
Print the elements that surround a givel geoel.
bool ComputeXInverse(TPZVec< REAL > &XD, TPZVec< REAL > &ksi, REAL Tol)
Computes the XInverse and returns true if ksi belongs to master element domain.
Definition: pzgeoel.cpp:686
Implements a geometric node in the pz environment. Geometry.
Definition: pzgnode.h:31
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...
int NSubElements() const
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
MElementType
Define the element types.
Definition: pzeltype.h:52
TPZGeoMesh & RefPatternMesh()
void Print(std::ostream &out=std::cout) const
Prints the useful information of a Refinement Pattern in a ostream file.
int Dimension() const
the dimension associated with the element/side
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.
string to_string(const string &value)
static void ModifyElementsBasedOnRefpFound(TPZAutoPointer< TPZRefPattern > &refpFound, TPZAutoPointer< TPZRefPattern > &refp, TPZVec< TPZGeoEl *> &elementVec)
Definition: pzeltype.h:61
int Side() const
Definition: pzgeoelside.h:169
static bool CompareTopologies(TPZAutoPointer< TPZRefPattern > refA, TPZAutoPointer< TPZRefPattern > refB, TPZTransform<> &fromAtoB, std::map< int, int > &pairNodes)
Returns if the given refPatterns (refA and refB) are topologicaly compatibles.
virtual TPZAutoPointer< TPZRefPattern > GetRefPattern() const
Returns the refinement pattern associated with the element.
Definition: pzgeoel.cpp:1716
static const std::string fNonInitializedName
Name for identifying a non initialized pattern.
Definition: TPZRefPattern.h:82
void SetMaterialId(int id)
Sets the material index of the element.
Definition: pzgeoel.h:395
void ComputeTransforms()
Calculates the transforms between the parametric coordinates of the sides of the son to the father&#39;s ...
static void TransformationTest(TPZRefPattern *refp)
Algorithm that evaluates the veracity of the hashings between sides of the elements children and corr...
virtual void SetRefPattern(TPZAutoPointer< TPZRefPattern >)
Defines the refinement pattern. It&#39;s used only in TPZGeoElRefPattern objects.
Definition: pzgeoel.cpp:1647
int Exists() const
Definition: pzgeoelside.h:175
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
static TPZAutoPointer< TPZRefPattern > ModelRefPattern(TPZGeoEl *gel, std::map< int, std::pair< TPZGeoEl *, std::map< int, int > > > &neighCorresp)
Returns the refpattern that matches the sides refinement by neighbours.
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
void GetCoordinates(TPZVec< REAL > &co)
Fill the coordinates of the node.
Definition: pzgnode.cpp:83
Definition: pzeltype.h:55
virtual void Print(std::ostream &out=std::cout)
Print all relevant data of the element to cout.
Definition: pzgeoel.cpp:238
void GenerateSideRefPatterns()
Generate the refinement patterns associated with the sides of the father element. ...
REAL ZeroTolerance()
Returns the tolerance to Zero value. Actually: .
Definition: pzreal.h:633
Contains the TPZVTKGeoMesh class which implements the graphical mesh to VTK environment to geometric ...
int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzfmatrix.h:548
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
int64_t SideNodeIndex(int nodenum) const
Returns the index of the nodenum node of side.
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
Contains the TPZRefPatternTools class which defines tools of pattern.
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
void ComputePartition()
It computers the partition of the sides of the father element using the sides of the children...
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138