NeoPZ
hdivCurvedJCompAppMath.cpp
Go to the documentation of this file.
1 //
2 // hdivCurvedJCompAppMath.cpp
3 // PZ
4 //
5 // Created by Douglas Castro on 18/06/15.
6 //
7 //
8 
9 #include "hdivCurvedJCompAppMath.h"
10 
11 const int norder = 6;
12 
13 bool hdivCurvedJCompAppMath::probAtCircle = false;
14 bool hdivCurvedJCompAppMath::probAtCylinder = false;
15 bool hdivCurvedJCompAppMath::probAtSphere = false;
16 
17 
18 hdivCurvedJCompAppMath::hdivCurvedJCompAppMath()
19 {
20 
21  fDim = 2;
22 
23  fmatId = 1;
24 
25  fdirichlet = 0;
26  fneumann = 1;
27 
28  fbc0 = -1;
29  fbc1 = -2;
30  fbc2 = -3;
31  fbc3 = -4;
32  fbc4 = -5;
33  fbc5 = -6;
34  fmatskeleton = -7;
35 
36  probAtCircle = false;
37  probAtCylinder = false;
38  probAtSphere = false;
39 
40  ftriang = false;
41 
42  isgeoblend = true;
43 }
44 
45 
46 hdivCurvedJCompAppMath::hdivCurvedJCompAppMath(geomDomain geodomain)
47 {
48  switch(geodomain){
49  case ECircle:
50  {
51  probAtCircle = true;
52  probAtCylinder = false;
53  probAtSphere = false;
54  }
55  break;
56  case ECylinder:
57  {
58  probAtCircle = false;
59  probAtCylinder = true;
60  probAtSphere = false;
61  }
62  break;
63  case ESphere:
64  {
65  probAtCircle = false;
66  probAtCylinder = false;
67  probAtSphere = true;
68  }
69  break;
70  default:
71  break;
72 
73  }
74 // switch(geodomain){
75 // case ECircle:
76 // {
77 // probAtCircle = 1;
78 // probAtCylinder = 0;
79 // probAtSphere = 0;
80 // }
81 // break;
82 // case ECylinder:
83 // {
84 // probAtCircle = 0;
85 // probAtCylinder = 1;
86 // probAtSphere = 0;
87 // }
88 // break;
89 // case ESphere:
90 // {
91 // probAtCircle = 0;
92 // probAtCylinder = 0;
93 // probAtSphere = 1;
94 // }
95 // break;
96 // default:
97 // break;
98 //
99 // }
100 }
101 
102 hdivCurvedJCompAppMath::~hdivCurvedJCompAppMath()
103 {
104 
105 }
106 
107 void hdivCurvedJCompAppMath::Run(geomDomain geodomain, ApproximationSpace problem, Eltype element, TPZVec<int> POrderBeginAndEnd, TPZVec<int> ndivinterval, TPZFMatrix< REAL > &errors)
108 {
109  if (POrderBeginAndEnd.size()!=2)
110  {
111  std::cout << " POrderBeginAndEnd Vector must have size == 2 and contain just two values, the first and last approximation orders. " << std::endl;
112  std::cout << " Some like POrderBeginAndEnd[0] = 1, POrderBeginAndEnd[1] = 3 " << std::endl;
113  return;
114  }
115 
116  if (ndivinterval.size()!=2)
117  {
118  std::cout << " ndivinterval Vector must have size == 2 and contain just two values, the first and last number of uniform refinements. " << std::endl;
119  std::cout << " Some like ndivinterval[0] = 0, ndivinterval[1] = 3 for simulations with 0, 1, 2 and 3 refinements " << std::endl;
120  std::cout << " Some like ndivinterval[0] = 1, ndivinterval[1] = 1 for simulations with 1 refinement " << std::endl;
121  return;
122  }
123 
124  int sizeErrorMatrix = (ndivinterval[1]-ndivinterval[0])*(POrderBeginAndEnd[1] - POrderBeginAndEnd[0]) + 1;
125  int errorposition = 0;
126 
127  errors.Resize(sizeErrorMatrix, 4);
128 
129  for(int p = POrderBeginAndEnd[0]; p <= POrderBeginAndEnd[1]; p++)
130  {
131 
132  for (int ndiv = ndivinterval[0]; ndiv<=ndivinterval[1]; ndiv++)
133  {
134  int ordemP = p;
135 
136  std::cout<< " BEGIN - Polinomial degree: " << ordemP << " rerfinement size (h): " << ndiv << std::endl;
137 
138  TPZGeoMesh *gmesh;
139  /*
140  switch (element) {
141  case EQuad:
142  {
143  gmesh = this->GMesh(fDim, false, ndiv);
144  }
145  break;
146  case ETriangle:
147  {
148  setTriangTrue();
149  gmesh = this->GMesh(fDim, true, ndiv);
150  }
151  break;
152  default:
153  break;
154  }*/
155 
156  hdivCurvedJCompAppMath geodomain;
157  if(probAtCircle)
158  {
159  gmesh = this->MakeCircle(ndiv);
160  }
161  else if(probAtCylinder)
162  {
163  gmesh = this->GMeshCilindricalMesh(ndiv);;
164 
165  }
166  else //if(probAtSphere)
167  {
168  gmesh = this->MakeSphereFromQuadrilateral(2, false, ndiv);
169  }
170 
171 // switch (geodomain) {
172 // case ECircle:
173 // {
174 // gmesh = this->MakeCircle(ndiv);
175 // }
176 // break;
177 // case ECylinder:
178 // {
179 // gmesh = this->GMeshCilindricalMesh(ndiv);;
180 // }
181 // break;
182 // case ESphere:
183 // {
184 // gmesh = this->MakeSphereFromQuadrilateral(2, false, ndiv);
185 // }
186 // break;
187 // default:
188 // break;
189 // }
190 
191  switch (problem) {
192  case EH1:
193  {
194  gmesh->SetDimension(fDim);
195  TPZCompMesh *cmeshH1 = this->CMeshH1(gmesh, ordemP, fDim);
196  //condensar
197  for (int64_t iel=0; iel<cmeshH1->NElements(); iel++) {
198  TPZCompEl *cel = cmeshH1->Element(iel);
199  if(!cel) continue;
200  new TPZCondensedCompEl(cel);
201  }
202 
203  cmeshH1->ExpandSolution();
204  cmeshH1->CleanUpUnconnectedNodes();
205 
206  TPZAnalysis anh1(cmeshH1, true);
207 
208  SolveSyst(anh1, cmeshH1);
209 
210  ErrorH1(cmeshH1, ordemP, ndiv, errorposition, errors);
211 
212  }
213  break;
214  case EHDiv:
215  {
216  DebugStop(); //Mudanca na topologia
217  }
218  break;
219  case EHDivStar:
220  {
221 
222  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
223 
224  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
225 
226  // P**
227  ChangeExternalOrderConnects(cmesh1);
228 
229 
230  //malha multifisica
231  TPZVec<TPZCompMesh *> meshvec(2);
232  meshvec[0] = cmesh1;
233  meshvec[1] = cmesh2;
234 
235  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
236 
237  TPZAnalysis an(mphysics, true);
238 
239  SolveSyst(an, mphysics);
240 
241  //Calculo do erro
242  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
243  TPZVec<REAL> erros;
244 
245  std::cout << "Computing Errors\n";
246 
247  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, errorposition, errors);
248  }
249  break;
250 
251  case EHDivStarStar:
252  {
253 
254  ordemP = ordemP + 1;
255 
256  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
257 
258  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
259  //P**
260  ChangeExternalOrderConnects(cmesh1);
261 
262  //malha multifisica
263  TPZVec<TPZCompMesh *> meshvec(2);
264  meshvec[0] = cmesh1;
265  meshvec[1] = cmesh2;
266 
267  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
268 
269  TPZAnalysis an(mphysics, true);
270 
271  SolveSyst(an, mphysics);
272 
273  //Calculo do erro
274  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
275  TPZVec<REAL> erros;
276 
277  std::cout << "Computing Errors\n";
278 
279  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, errorposition, errors);
280  }
281  break;
282 
283  default:
284  break;
285  }
286 
287  std::cout<< " END - polinomial degree: " << ordemP << " refinement size (h): " << ndiv << std::endl;
288 
289  errorposition++;
290 
291  }
292  }
293 
294  std::cout<< " The End " << std::endl;
295 
296 
297 }
298 
299 void hdivCurvedJCompAppMath::PrintErrors(geomDomain geodomain, ApproximationSpace problem, Eltype element, TPZVec<int> POrderBeginAndEnd, TPZVec<int> ndivinterval, TPZVec<REAL> &errors, std::ostream &output)
300 {
301  if (POrderBeginAndEnd.size()!=2)
302  {
303  std::cout << " POrderBeginAndEnd Vector must have size == 2 and contain just two values, the first and last approximation orders. " << std::endl;
304  std::cout << " Some like POrderBeginAndEnd[0] = 1, POrderBeginAndEnd[1] = 3 " << std::endl;
305  return;
306  }
307 
308  if (ndivinterval.size()!=2)
309  {
310  std::cout << " ndivinterval Vector must have size == 2 and contain just two values, the first and last number of uniform refinements. " << std::endl;
311  std::cout << " Some like ndivinterval[0] = 0, ndivinterval[1] = 3 for simulations with 0, 1, 2 and 3 refinements " << std::endl;
312  std::cout << " Some like ndivinterval[0] = 1, ndivinterval[1] = 1 for simulations with 1 refinement " << std::endl;
313  return;
314  }
315 
316  for(int p = POrderBeginAndEnd[0]; p <= POrderBeginAndEnd[1]; p++)
317  {
318  output << "\n WHEN p = " << p << " \n " << endl;
319  output << "ndiv " << setw(6) << "DoFTot" << setw(20) << "DofCond" << setw(28) << "PrimalL2Error" << setw(35) << "L2DualError" << endl;
320 
321  for (int ndiv = ndivinterval[0]; ndiv<=ndivinterval[1]; ndiv++)
322  {
323  int ordemP = p;
324 
325  std::cout<< " BEGIN (Quadrilateral Domain) - Polinomial degree: " << ordemP << " rerfinement size (h): " << ndiv << std::endl;
326 
327  TPZGeoMesh *gmesh;
328 
329  hdivCurvedJCompAppMath geodomain;
330 
331  if(probAtCircle)
332  {
333  gmesh = this->MakeCircle(ndiv);
334  }
335  else if(probAtCylinder)
336  {
337  gmesh = this->GMeshCilindricalMesh(ndiv);;
338 
339  }
340  else //if(probAtSphere)
341  {
342  gmesh = this->MakeSphereFromQuadrilateral(2, false, ndiv);
343  }
344 // switch (geodomain) {
345 // case ECircle:
346 // {
347 // gmesh = this->MakeCircle( ndiv);
348 // }
349 // break;
350 // case ECylinder:
351 // {
352 // gmesh = this->GMeshCilindricalMesh(ndiv);;
353 // }
354 // break;
355 // case ESphere:
356 // {
357 // gmesh = this->MakeSphereFromQuadrilateral(2, false, ndiv);
358 // }
359 // break;
360 // default:
361 // break;
362 // }
363 
364  gmesh->SetDimension(fDim);
365 
366  switch (problem) {
367  case EH1:
368  {
369 
370  TPZCompMesh *cmeshH1 = this->CMeshH1(gmesh, ordemP, fDim);
371 
372  int dofTotal = cmeshH1->NEquations();
373 
374  //condensar
375  for (int64_t iel=0; iel<cmeshH1->NElements(); iel++) {
376  TPZCompEl *cel = cmeshH1->Element(iel);
377  if(!cel) continue;
378  new TPZCondensedCompEl(cel);
379  }
380 
381  cmeshH1->ExpandSolution();
382  cmeshH1->CleanUpUnconnectedNodes();
383 
384  int dofCondensed = cmeshH1->NEquations();
385 
386  TPZAnalysis anh1(cmeshH1, true);
387 
388  SolveSyst(anh1, cmeshH1);
389 
390  ErrorH1(cmeshH1, ordemP, ndiv, output, dofTotal, dofCondensed);
391 
392  }
393  break;
394  case EHDiv:
395  {
396  DebugStop(); //Mudanca na topologia
397  }
398  break;
399  case EHDivStar:
400  {
401 
402  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
403 
404  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
405 
406  // para rodar P**
407  ChangeExternalOrderConnects(cmesh1);
408 
409  int DofCond, DoFT;
410  DoFT = cmesh1->NEquations() + cmesh2->NEquations();
411 
412  //malha multifisica
413  TPZVec<TPZCompMesh *> meshvec(2);
414  meshvec[0] = cmesh1;
415  meshvec[1] = cmesh2;
416 
417  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
418 
419  DofCond = mphysics->NEquations();
420 
421  TPZAnalysis an(mphysics, true);
422 
423  SolveSyst(an, mphysics);
424 
425  //Calculo do erro
426  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
427  TPZVec<REAL> erros;
428 
429  std::cout << "Computing Errors\n";
430 
431  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, output, DoFT, DofCond);
432  }
433  break;
434 
435  case EHDivStarStar:
436  {
437 
438 
439  ordemP = ordemP + 1;
440 
441  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
442 
443  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
444 
445  // para rodar P**
446  ChangeExternalOrderConnects(cmesh1);
447 
448  int DofCond, DoFT;
449  DoFT = cmesh1->NEquations() + cmesh2->NEquations();
450 
451  //malha multifisica
452  TPZVec<TPZCompMesh *> meshvec(2);
453  meshvec[0] = cmesh1;
454  meshvec[1] = cmesh2;
455 
456  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
457 
458  DofCond = mphysics->NEquations();
459 
460  TPZAnalysis an(mphysics, true);
461 
462  SolveSyst(an, mphysics);
463 
464  //Calculo do erro
465  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
466  TPZVec<REAL> erros;
467 
468  std::cout << "Computing Errors\n";
469 
470  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, output, DoFT, DofCond);
471  }
472  break;
473 
474  default:
475  break;
476  }
477 
478 
479  }
480 
481  output << "\n ----------------------------------------------------------------------------- " << endl;
482  std::cout<< " END (Quadrilateral Domain) - Polinomial degree: " << p << std::endl;
483 
484  }
485 
486  std::cout<< " The End " << std::endl;
487 
488 
489 }
490 
491 TPZGeoMesh *hdivCurvedJCompAppMath::MakeCircle( int ndiv)
492 {
493 
494  if( fDim!= 2)
495  {
496  DebugStop();
497  }
498 
499  TPZGeoMesh * geomesh = new TPZGeoMesh;
500 
501  int nodes = 17 + 8;
502  REAL radius = 1.0;
503  REAL innerradius = radius/2.0;
504  geomesh->SetMaxNodeId(nodes-1);
505  geomesh->NodeVec().Resize(nodes);
506  TPZManVector<TPZGeoNode,7> Node(nodes);
507 
508  TPZManVector<int64_t,8> TopolQQuadrilateral(8);
509  TPZManVector<int64_t,8> TopolQuadrilateral(4);
510  TPZManVector<int64_t,6> TopolQTriangle(6);
511  TPZManVector<int64_t,2> TopolLine(2);
512  TPZManVector<int64_t,3> TopolArc(3);
513  TPZManVector<REAL,3> coord(3,0.);
514  TPZVec<REAL> xc(3,0.);
515 
516 
517  int nodeindex = 0;
518 
519  for (int inode = 0; inode < 8 ; inode++) {
520  // i node
521  coord = ParametricCircle(innerradius, inode * M_PI/4.0);
522  geomesh->NodeVec()[nodeindex].SetCoord(coord);
523  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
524  nodeindex++;
525  }
526 
527  for (int inode = 0; inode < 8 ; inode++) {
528  // i node
529  coord = ParametricCircle(radius, inode * M_PI/4.0);
530  geomesh->NodeVec()[nodeindex].SetCoord(coord);
531  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
532  nodeindex++;
533  }
534 
535  for (int inode = 0; inode < 4 ; inode++) {
536  // i node
537  coord = ParametricCircle(0.5*innerradius, inode * M_PI/2.0);
538  geomesh->NodeVec()[nodeindex].SetCoord(coord);
539  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
540  nodeindex++;
541  }
542 
543  for (int inode = 0; inode < 4 ; inode++) {
544  // i node
545  coord = ParametricCircle(1.5*innerradius, inode * M_PI/2.0);
546  geomesh->NodeVec()[nodeindex].SetCoord(coord);
547  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
548  nodeindex++;
549  }
550 
551  // 24 node id at circle center
552  coord = ParametricCircle(0.0,0.0);
553  geomesh->NodeVec()[nodeindex].SetCoord(coord);
554  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
555  nodeindex++;
556 
557  int elementid = 0;
558 
559  // outer domain
560 
561  TopolQuadrilateral[0] = 0;
562  TopolQuadrilateral[1] = 8;
563  TopolQuadrilateral[2] = 10;
564  TopolQuadrilateral[3] = 2;
565  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend<pzgeom::TPZGeoQuad > > (elementid,TopolQuadrilateral, fmatId,*geomesh);
566  elementid++;
567 
568  TopolQuadrilateral[0] = 2;
569  TopolQuadrilateral[1] = 10;
570  TopolQuadrilateral[2] = 12;
571  TopolQuadrilateral[3] = 4;
572  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend<pzgeom::TPZGeoQuad > > (elementid,TopolQuadrilateral, fmatId,*geomesh);
573  elementid++;
574 
575  TopolQuadrilateral[0] = 4;
576  TopolQuadrilateral[1] = 12;
577  TopolQuadrilateral[2] = 14;
578  TopolQuadrilateral[3] = 6;
579 
580  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend<pzgeom::TPZGeoQuad > > (elementid,TopolQuadrilateral, fmatId,*geomesh);
581  elementid++;
582 
583  TopolQuadrilateral[0] = 6;
584  TopolQuadrilateral[1] = 14;
585  TopolQuadrilateral[2] = 8;
586  TopolQuadrilateral[3] = 0;
587  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend<pzgeom::TPZGeoQuad > > (elementid,TopolQuadrilateral, fmatId,*geomesh);
588  elementid++;
589 
590 
591  // inner domain
592 
593 
594  TopolQTriangle[0] = 0;
595  TopolQTriangle[1] = 2;
596  TopolQTriangle[2] = 24;
597  TopolQTriangle[3] = 1;
598  TopolQTriangle[4] = 17;
599  TopolQTriangle[5] = 16;
600  new TPZGeoElRefPattern< pzgeom::TPZQuadraticTrig > (elementid,TopolQTriangle, fmatId,*geomesh);
601  elementid++;
602 
603  TopolQTriangle[0] = 2;
604  TopolQTriangle[1] = 4;
605  TopolQTriangle[2] = 24;
606  TopolQTriangle[3] = 3;
607  TopolQTriangle[4] = 18;
608  TopolQTriangle[5] = 17;
609  new TPZGeoElRefPattern< pzgeom::TPZQuadraticTrig > (elementid,TopolQTriangle, fmatId,*geomesh);
610  elementid++;
611 
612  TopolQTriangle[0] = 4;
613  TopolQTriangle[1] = 6;
614  TopolQTriangle[2] = 24;
615  TopolQTriangle[3] = 5;
616  TopolQTriangle[4] = 19;
617  TopolQTriangle[5] = 18;
618  new TPZGeoElRefPattern< pzgeom::TPZQuadraticTrig > (elementid,TopolQTriangle, fmatId,*geomesh);
619  elementid++;
620 
621  TopolQTriangle[0] = 6;
622  TopolQTriangle[1] = 0;
623  TopolQTriangle[2] = 24;
624  TopolQTriangle[3] = 7;
625  TopolQTriangle[4] = 16;
626  TopolQTriangle[5] = 19;
627  new TPZGeoElRefPattern< pzgeom::TPZQuadraticTrig > (elementid,TopolQTriangle, fmatId,*geomesh);
628  elementid++;
629 
630 
631 
632  // outer arcs bc's
633 
634  TopolArc[0] = 8;
635  TopolArc[1] = 10;
636  TopolArc[2] = 9;
637  new TPZGeoElRefPattern< pzgeom::TPZArc3D > (elementid,TopolArc, fbc1,*geomesh);
638  elementid++;
639 
640  TopolArc[0] = 10;
641  TopolArc[1] = 12;
642  TopolArc[2] = 11;
643  new TPZGeoElRefPattern< pzgeom::TPZArc3D > (elementid,TopolArc, fbc1,*geomesh);
644  elementid++;
645 
646  TopolArc[0] = 12;
647  TopolArc[1] = 14;
648  TopolArc[2] = 13;
649  new TPZGeoElRefPattern< pzgeom::TPZArc3D > (elementid,TopolArc, fbc1,*geomesh);
650  elementid++;
651 
652  TopolArc[0] = 14;
653  TopolArc[1] = 8;
654  TopolArc[2] = 15;
655  new TPZGeoElRefPattern< pzgeom::TPZArc3D > (elementid,TopolArc, fbc1,*geomesh);
656  elementid++;
657 
658 
659 
660  geomesh->BuildConnectivity();
661 
662  int nref = ndiv;
663  TPZVec<TPZGeoEl *> sons;
664  for (int iref = 0; iref < nref; iref++) {
665  int nel = geomesh->NElements();
666  for (int iel = 0; iel < nel; iel++) {
667  TPZGeoEl *gel = geomesh->ElementVec()[iel];
668  if (gel->HasSubElement()) {
669  continue;
670  }
671  gel->Divide(sons);
672  }
673  }
674 
675 
676  std::ofstream out("CircleMixed.vtk");
677  TPZVTKGeoMesh::PrintGMeshVTK(geomesh, out, true);
678 
679  return geomesh;
680 }
681 
682 TPZManVector<REAL,3> hdivCurvedJCompAppMath::ParametricCircle(REAL radius,REAL theta)
683 {
684  TPZManVector<REAL,3> xcoor(3,0.0);
685  xcoor[0] = radius * cos(theta);
686  xcoor[1] = radius * sin(theta);
687  xcoor[2] = 0.0 ;
688  return xcoor;
689 }
690 
691 TPZGeoMesh *hdivCurvedJCompAppMath::MakeSphereFromQuadrilateral(int dimensao, bool triang, int ndiv)
692 {
693  TPZGeoMesh * geomesh = new TPZGeoMesh;
694 
695  int nodes = 8;
696  REAL radius = 1.0;
697  geomesh->SetMaxNodeId(nodes-1);
698  geomesh->NodeVec().Resize(nodes);
699  TPZManVector<TPZGeoNode,4> Node(nodes);
700 
701  TPZManVector<int64_t,2> TopolQuad(4);
702  TPZManVector<int64_t,1> TopolPoint(1);
703  TPZManVector<int64_t,2> TopolLine(2);
704  TPZManVector<REAL,3> coord(3,0.);
705  TPZVec<REAL> xc(3,0.);
706 
707  REAL cphi = atan(sqrt(2.0));
708 
709  int nodeindex = 0;
710 
711  coord = ParametricSphere(radius,M_PI-cphi,M_PI/4.0);
712  geomesh->NodeVec()[nodeindex].SetCoord(coord);
713  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
714  nodeindex++;
715 
716  coord = ParametricSphere(radius,cphi,M_PI/4.0);
717  geomesh->NodeVec()[nodeindex].SetCoord(coord);
718  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
719  nodeindex++;
720 
721  coord = ParametricSphere(radius,cphi,-M_PI/4.0);
722  geomesh->NodeVec()[nodeindex].SetCoord(coord);
723  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
724  nodeindex++;
725 
726  coord = ParametricSphere(radius,M_PI-cphi,-M_PI/4.0);
727  geomesh->NodeVec()[nodeindex].SetCoord(coord);
728  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
729  nodeindex++;
730 
731 
732  coord = ParametricSphere(radius,M_PI-cphi,3.0*M_PI/4.0);
733  geomesh->NodeVec()[nodeindex].SetCoord(coord);
734  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
735  nodeindex++;
736 
737  coord = ParametricSphere(radius,cphi,3.0*M_PI/4.0);
738  geomesh->NodeVec()[nodeindex].SetCoord(coord);
739  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
740  nodeindex++;
741 
742  coord = ParametricSphere(radius,cphi,-3.0*M_PI/4.0);
743  geomesh->NodeVec()[nodeindex].SetCoord(coord);
744  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
745  nodeindex++;
746 
747  coord = ParametricSphere(radius,M_PI-cphi,-3.0*M_PI/4.0);
748  geomesh->NodeVec()[nodeindex].SetCoord(coord);
749  geomesh->NodeVec()[nodeindex].SetNodeId(nodeindex);
750  nodeindex++;
751 
752 
753  int id = 0;
754  int matid = 1;
755 
756  TopolQuad[0] = 0;
757  TopolQuad[1] = 1;
758  TopolQuad[2] = 2;
759  TopolQuad[3] = 3;
760  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad1 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
761  quad1->Geom().SetData(radius, xc);
762  id++;
763 
764  TopolQuad[0] = 4;
765  TopolQuad[1] = 5;
766  TopolQuad[2] = 6;
767  TopolQuad[3] = 7;
768  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad2 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
769  quad2->Geom().SetData(radius, xc);
770  id++;
771 
772  TopolQuad[0] = 0;
773  TopolQuad[1] = 4;
774  TopolQuad[2] = 5;
775  TopolQuad[3] = 1;
776  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad3 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
777  quad3->Geom().SetData(radius, xc);
778  id++;
779 
780  TopolQuad[0] = 3;
781  TopolQuad[1] = 7;
782  TopolQuad[2] = 6;
783  TopolQuad[3] = 2;
784  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad4 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
785  quad4->Geom().SetData(radius, xc);
786  id++;
787 
788  TopolQuad[0] = 0;
789  TopolQuad[1] = 4;
790  TopolQuad[2] = 7;
791  TopolQuad[3] = 3;
792  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad5 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
793  quad5->Geom().SetData(radius, xc);
794  id++;
795 
796  TopolQuad[0] = 1;
797  TopolQuad[1] = 5;
798  TopolQuad[2] = 6;
799  TopolQuad[3] = 2;
800  TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > * quad6 = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoQuad > > (id,TopolQuad,matid,*geomesh);
801  quad6->Geom().SetData(radius, xc);
802  id++;
803 
804  TopolLine[0] = 0;
805  TopolLine[1] = 4;
806  new TPZGeoElRefPattern < pzgeom::TPZGeoBlend<pzgeom::TPZGeoLinear> > (id,TopolLine, fbc0, *geomesh);
807  id++;
808 
809  TopolLine[0] = 0;
810  TopolLine[1] = 4;
811  new TPZGeoElRefPattern < pzgeom::TPZGeoBlend<pzgeom::TPZGeoLinear> > (id,TopolLine, fbc1, *geomesh);
812  id++;
813 
814 
815 // TopolLine[0] = 0;
816 // TopolLine[1] = 4;
817 // TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoLinear> > * arcleft = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoLinear> > (id,TopolLine, fbc0, *geomesh);
818 // arcleft->Geom().SetData(radius, xc);
819 // id++;
820 //
821 // TopolLine[0] = 0;
822 // TopolLine[1] = 4;
823 // TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoLinear> > * arcrigth = new TPZGeoElRefPattern< pzgeom::TPZQuadSphere< pzgeom::TPZGeoLinear> > (id,TopolLine, fbc1, *geomesh);
824 // arcrigth->Geom().SetData(radius, xc);
825 // id++;
826 
827  geomesh->BuildConnectivity();
828 
829  TPZVec<TPZGeoEl *> sons;
830  const int nref = ndiv;
831  for (int iref = 0; iref < nref; iref++) {
832  int nel = geomesh->NElements();
833  for (int iel = 0; iel < nel; iel++) {
834  TPZGeoEl *gel = geomesh->ElementVec()[iel];
835  if(!gel->HasSubElement())
836  {
837  gel->Divide(sons);
838  }
839  }
840  }
841 
842  int axis = 1;
843  REAL angle = -45.0;
844  this->RotateGeomesh(geomesh, angle, axis);
845 
846 // TPZCheckGeom * checkgeo = new TPZCheckGeom(geomesh);
847 // int badmeshQ = checkgeo->PerformCheck();
848 //
849 // if (badmeshQ) {
850 // DebugStop();
851 // }
852 
853 // std::ofstream outfile("NiceSphere.vtk");
854 // TPZVTKGeoMesh::PrintGMeshVTK(geomesh, outfile, true);
855 //
856  return geomesh;
857 }
858 
859 void hdivCurvedJCompAppMath::RotateGeomesh(TPZGeoMesh *gmesh, REAL CounterClockwiseAngle, int &Axis)
860 {
861  REAL theta = (M_PI/180.0)*CounterClockwiseAngle;
862  // It represents a 3D rotation around the z axis.
863  TPZFMatrix<STATE> RotationMatrix(3,3,0.0);
864 
865  switch (Axis) {
866  case 1:
867  {
868  RotationMatrix(0,0) = 1.0;
869  RotationMatrix(1,1) = +cos(theta);
870  RotationMatrix(1,2) = -sin(theta);
871  RotationMatrix(2,1) = +sin(theta);
872  RotationMatrix(2,2) = +cos(theta);
873  }
874  break;
875  case 2:
876  {
877  RotationMatrix(0,0) = +cos(theta);
878  RotationMatrix(0,2) = +sin(theta);
879  RotationMatrix(1,1) = 1.0;
880  RotationMatrix(2,0) = -sin(theta);
881  RotationMatrix(2,2) = +cos(theta);
882  }
883  break;
884  case 3:
885  {
886  RotationMatrix(0,0) = +cos(theta);
887  RotationMatrix(0,1) = -sin(theta);
888  RotationMatrix(1,0) = +sin(theta);
889  RotationMatrix(1,1) = +cos(theta);
890  RotationMatrix(2,2) = 1.0;
891  }
892  break;
893  default:
894  {
895  RotationMatrix(0,0) = +cos(theta);
896  RotationMatrix(0,1) = -sin(theta);
897  RotationMatrix(1,0) = +sin(theta);
898  RotationMatrix(1,1) = +cos(theta);
899  RotationMatrix(2,2) = 1.0;
900  }
901  break;
902  }
903 
904  TPZVec<REAL> iCoords(3,0.0);
905  TPZVec<REAL> iCoordsRotated(3,0.0);
906 
907  //RotationMatrix.Print("Rotation = ");
908 
909  int NumberofGeoNodes = gmesh->NNodes();
910  for (int inode = 0; inode < NumberofGeoNodes; inode++)
911  {
912  TPZGeoNode GeoNode = gmesh->NodeVec()[inode];
913  GeoNode.GetCoordinates(iCoords);
914  // Apply rotation
915  iCoordsRotated[0] = RotationMatrix(0,0)*iCoords[0]+RotationMatrix(0,1)*iCoords[1]+RotationMatrix(0,2)*iCoords[2];
916  iCoordsRotated[1] = RotationMatrix(1,0)*iCoords[0]+RotationMatrix(1,1)*iCoords[1]+RotationMatrix(1,2)*iCoords[2];
917  iCoordsRotated[2] = RotationMatrix(2,0)*iCoords[0]+RotationMatrix(2,1)*iCoords[1]+RotationMatrix(2,2)*iCoords[2];
918  GeoNode.SetCoord(iCoordsRotated);
919  gmesh->NodeVec()[inode] = GeoNode;
920  }
921 
922 }
923 
924 
925 
926 TPZManVector<REAL,3> hdivCurvedJCompAppMath::ParametricSphere(REAL radius,REAL phi,REAL theta)
927 {
928  TPZManVector<REAL,3> xcoor(3,0.0);
929  xcoor[0] = radius * cos(theta) * sin(phi);
930  xcoor[1] = radius * sin(theta) * sin(phi);
931  xcoor[2] = radius * cos(phi) ;
932  return xcoor;
933 }
934 
935 TPZGeoMesh *hdivCurvedJCompAppMath::GMeshCilindricalMesh( int ndiv)
936 {
938  TPZGeoMesh * gmesh = new TPZGeoMesh;
939 
941 
942  int nodenumber = 10;
943  REAL r = 1.0;
944  // o angulo theta varia de theta0 a theta1
945  REAL theta0 = 0.0;
946  REAL theta1 = M_PI;
947  // z varia de z0 a z1
948  REAL z0 = -M_PI/2.0;
949  REAL z1 = M_PI/2.0;
950 
951  TPZVec<REAL> coord(3,0.0);
952  int Axis = 3;
953  REAL angrot = 0.0;
954 
955  gmesh = new TPZGeoMesh;
956  gmesh->NodeVec().Resize(nodenumber);
957 
958  // Setting node coordinates
959  int id = 0;
960  //0
961  coord[0] = r*cos(theta0);
962  coord[1] = r*sin(theta0);
963  coord[2] = z0;
964  RotateNode(coord, angrot, Axis);
965 
966  gmesh->NodeVec()[id].SetNodeId(id);
967  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
968  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
969  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
970  id++;
971 
972  //1
973  coord[0] = r*cos(theta1);
974  coord[1] = r*sin(theta1);
975  coord[2] = z0;
976  RotateNode(coord, angrot, Axis);
977 
978  gmesh->NodeVec()[id].SetNodeId(id);
979  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
980  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
981  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
982  id++;
983 
984  //2
985  coord[0] = r*cos(theta1);
986  coord[1] = r*sin(theta1);
987  coord[2] = z1;
988  RotateNode(coord, angrot, Axis);
989 
990  gmesh->NodeVec()[id].SetNodeId(id);
991  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
992  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
993  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
994  id++;
995 
996  //3
997  coord[0] = r*cos(theta0);
998  coord[1] = r*sin(theta0);
999  coord[2] = z1;
1000  RotateNode(coord, angrot, Axis);
1001 
1002  gmesh->NodeVec()[id].SetNodeId(id);
1003  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1004  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1005  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1006  id++;
1007 
1008  //4
1009  coord[0] = r*cos(M_PI/4.0);
1010  coord[1] = r*sin(M_PI/4.0);
1011  coord[2] = z0;
1012  RotateNode(coord, angrot, Axis);
1013 
1014  gmesh->NodeVec()[id].SetNodeId(id);
1015  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1016  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1017  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1018  id++;
1019 
1020  //5
1021  coord[0] = r*cos(M_PI/2.0);
1022  coord[1] = r*sin(M_PI/2.0);
1023  coord[2] = z0;
1024  RotateNode(coord, angrot, Axis);
1025 
1026  gmesh->NodeVec()[id].SetNodeId(id);
1027  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1028  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1029  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1030  id++;
1031 
1032  //6
1033  coord[0] = r*cos(3.0*M_PI/4.0);
1034  coord[1] = r*sin(3.0*M_PI/4.0);
1035  coord[2] = z0;
1036  RotateNode(coord, angrot, Axis);
1037 
1038  gmesh->NodeVec()[id].SetNodeId(id);
1039  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1040  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1041  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1042  id++;
1043 
1044  //7
1045  coord[0] = r*cos(M_PI/4.0);
1046  coord[1] = r*sin(M_PI/4.0);
1047  coord[2] = z1;
1048  RotateNode(coord, angrot, Axis);
1049 
1050  gmesh->NodeVec()[id].SetNodeId(id);
1051  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1052  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1053  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1054  id++;
1055 
1056  //8
1057  coord[0] = r*cos(M_PI/2.0);
1058  coord[1] = r*sin(M_PI/2.0);
1059  coord[2] = z1;
1060  RotateNode(coord, angrot, Axis);
1061 
1062  gmesh->NodeVec()[id].SetNodeId(id);
1063  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1064  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1065  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1066  id++;
1067 
1068  //9
1069  coord[0] = r*cos(3.0*M_PI/4.0);
1070  coord[1] = r*sin(3.0*M_PI/4.0);
1071  coord[2] = z1;
1072  RotateNode(coord, angrot, Axis);
1073 
1074  gmesh->NodeVec()[id].SetNodeId(id);
1075  gmesh->NodeVec()[id].SetCoord(0,coord[0]);//coord X
1076  gmesh->NodeVec()[id].SetCoord(1,coord[1]);//coord Y
1077  gmesh->NodeVec()[id].SetCoord(2,coord[2]);//coord Z
1078  id++;
1079 
1080  int elementid = 0;
1081  TPZVec < int64_t > nodeindex(3,0L);
1082 
1083  // Definition of Arc coordenates
1084 
1085 
1086  nodeindex.resize(2);
1087 
1088  // Create Geometrical Arc #2
1089  nodeindex[0] = 0;
1090  nodeindex[1] = 3;
1091  new TPZGeoElRefPattern < pzgeom::TPZGeoLinear > (elementid,nodeindex, fbc1, *gmesh);
1092  elementid++;
1093 
1094  nodeindex.resize(3);
1095  // Create Geometrical Arc #1
1096  nodeindex[0] = 0;
1097  nodeindex[1] = 5;
1098  nodeindex[2] = 4;
1099  new TPZGeoElRefPattern < pzgeom::TPZArc3D > (elementid,nodeindex, fbc2, *gmesh);
1100  elementid++;
1101  // Create Geometrical Arc #1
1102  nodeindex[0] = 5;
1103  nodeindex[1] = 1;
1104  nodeindex[2] = 6;
1105  new TPZGeoElRefPattern < pzgeom::TPZArc3D > (elementid,nodeindex, fbc2, *gmesh);
1106  elementid++;
1107 
1108 
1109  nodeindex.resize(2);
1110  // Create Geometrical Arc #4
1111  nodeindex[0] = 1;
1112  nodeindex[1] = 2;
1113  new TPZGeoElRefPattern < pzgeom::TPZGeoLinear > (elementid,nodeindex, fbc3, *gmesh);
1114 
1115 
1116  nodeindex.resize(3);
1117 
1118  // Create Geometrical Arc #3
1119  nodeindex[0] = 3;
1120  nodeindex[1] = 8;
1121  nodeindex[2] = 7;
1122  new TPZGeoElRefPattern < pzgeom::TPZArc3D > (elementid,nodeindex, fbc4, *gmesh);
1123  elementid++;
1124  // Create Geometrical Arc #1
1125  nodeindex[0] = 8;
1126  nodeindex[1] = 2;
1127  nodeindex[2] = 9;
1128  new TPZGeoElRefPattern < pzgeom::TPZArc3D > (elementid,nodeindex, fbc2, *gmesh);
1129  elementid++;
1130 
1131 
1132  nodeindex.resize(4);
1133 
1134  // Create Geometrical Quad #1
1135  nodeindex[0] = 0;
1136  nodeindex[1] = 5;
1137  nodeindex[2] = 8;
1138  nodeindex[3] = 3;
1139  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend < pzgeom::TPZGeoQuad > > (elementid,nodeindex, fmatId,*gmesh);
1140  elementid++;
1141  // Create Geometrical Quad #1
1142  nodeindex[0] = 1;
1143  nodeindex[1] = 2;
1144  nodeindex[2] = 8;
1145  nodeindex[3] = 5;
1146  new TPZGeoElRefPattern< pzgeom::TPZGeoBlend < pzgeom::TPZGeoQuad > > (elementid,nodeindex, fmatId,*gmesh);
1147  elementid++;
1148 
1149 
1150 
1151 
1152  //CONCLUINDO A CONSTRUCAO DA MALHA GEOMETRICA
1153  gmesh->BuildConnectivity();
1154 
1155  TPZVec<TPZGeoEl *> sons;
1156  const int nref = ndiv;
1157  for (int iref = 0; iref < nref; iref++) {
1158  int nel = gmesh->NElements();
1159  for (int iel = 0; iel < nel; iel++) {
1160  TPZGeoEl *gel = gmesh->ElementVec()[iel];
1161  if(!gel->HasSubElement()) gel->Divide(sons);
1162  }
1163  }
1164 
1165 // std::ofstream outfile("malhaCilindrica.vtk");
1166 // TPZVTKGeoMesh::PrintGMeshVTK(gmesh, outfile, true);
1167 
1168  return gmesh;
1169 
1170 }
1171 
1172 void hdivCurvedJCompAppMath::RotateNode(TPZVec<REAL> &iCoords, REAL CounterClockwiseAngle, int &Axis)
1173 {
1174  REAL theta = (M_PI/180.0)*CounterClockwiseAngle;
1175  // It represents a 3D rotation around the z axis.
1176  TPZFMatrix<REAL> RotationMatrix(3,3,0.0);
1177 
1178  switch (Axis) {
1179  case 1:
1180  {
1181  RotationMatrix(0,0) = 1.0;
1182  RotationMatrix(1,1) = +cos(theta);
1183  RotationMatrix(1,2) = -sin(theta);
1184  RotationMatrix(2,1) = +sin(theta);
1185  RotationMatrix(2,2) = +cos(theta);
1186  }
1187  break;
1188  case 2:
1189  {
1190  RotationMatrix(0,0) = +cos(theta);
1191  RotationMatrix(0,2) = +sin(theta);
1192  RotationMatrix(1,1) = 1.0;
1193  RotationMatrix(2,0) = -sin(theta);
1194  RotationMatrix(2,2) = +cos(theta);
1195  }
1196  break;
1197  case 3:
1198  {
1199  RotationMatrix(0,0) = +cos(theta);
1200  RotationMatrix(0,1) = -sin(theta);
1201  RotationMatrix(1,0) = +sin(theta);
1202  RotationMatrix(1,1) = +cos(theta);
1203  RotationMatrix(2,2) = 1.0;
1204  }
1205  break;
1206  default:
1207  {
1208  RotationMatrix(0,0) = +cos(theta);
1209  RotationMatrix(0,1) = -sin(theta);
1210  RotationMatrix(1,0) = +sin(theta);
1211  RotationMatrix(1,1) = +cos(theta);
1212  RotationMatrix(2,2) = 1.0;
1213  }
1214  break;
1215  }
1216 
1217  TPZVec<REAL> iCoordsRotated(3,0.0);
1218  // Apply rotation
1219  iCoordsRotated[0] = RotationMatrix(0,0)*iCoords[0]+RotationMatrix(0,1)*iCoords[1]+RotationMatrix(0,2)*iCoords[2];
1220  iCoordsRotated[1] = RotationMatrix(1,0)*iCoords[0]+RotationMatrix(1,1)*iCoords[1]+RotationMatrix(1,2)*iCoords[2];
1221  iCoordsRotated[2] = RotationMatrix(2,0)*iCoords[0]+RotationMatrix(2,1)*iCoords[1]+RotationMatrix(2,2)*iCoords[2];
1222  iCoords = iCoordsRotated;
1223 }
1224 
1225 
1226 void hdivCurvedJCompAppMath::SolExata(const TPZVec<REAL> &pt, TPZVec<STATE> &solp, TPZFMatrix<STATE> &flux){
1227 
1228  solp.resize(1);
1229  solp[0]=0.;
1230  flux.Resize(3,1);
1231 
1232  solp[0]=0.;
1233  flux(0,0)=0.0;
1234  flux(1,0)=0.0;
1235  flux(2,0)=0.0;
1236 
1237  REAL x = pt[0];
1238  REAL y = pt[1];
1239  REAL z = pt[2];
1240 
1241  if(probAtCircle)
1242  {
1243  //+======================================================================================
1244  // Circle
1245 
1246  REAL r = sqrt(x*x+y*y+z*z);
1247  REAL theta = atan2(y,x);
1248 
1249  REAL costheta = cos(theta);
1250  REAL sintheta = sin(theta);
1251  REAL cos2theta = cos(2.0*theta);
1252  REAL sin2theta = sin(2.0*theta);
1253  REAL r4 = r*r*r*r;
1254 
1255 
1256 
1257  solp[0] = r4*cos2theta*sin2theta;
1258 
1259  // Gradient computations
1260 
1261  REAL Runitx = costheta;
1262  REAL Runity = sintheta;
1263  REAL Runitz = 0.0;
1264 
1265  REAL Thetaunitx = -sintheta;
1266  REAL Thetaunity = costheta;
1267  REAL Thetaunitz = 0.0;
1268 
1269  REAL Zunitx = 0.0;
1270  REAL Zunity = 0.0;
1271  REAL Zunitz = 1.0;
1272 
1273  REAL dfdR = 4.0*r*r*r*cos2theta*sin2theta;
1274  REAL dfdTheta = (1.0)*(2.0*r*r*r*(cos2theta*cos2theta-sin2theta*sin2theta));
1275  REAL dfdZ = 0.0;
1276 
1277  flux(0,0) = -1.0*(dfdR * Runitx + dfdTheta * Thetaunitx + dfdZ * Zunitx);
1278  flux(1,0) = -1.0*(dfdR * Runity + dfdTheta * Thetaunity + dfdZ * Zunity);
1279  flux(2,0) = -1.0*(dfdR * Runitz + dfdTheta * Thetaunitz + dfdZ * Zunitz);
1280  }
1281  else if(probAtCylinder)
1282  {
1283  //+======================================================================================
1284  // Cylinder
1285  REAL r = sqrt(x*x+y*y);
1286  REAL theta = atan2(y,x);
1287 
1288  solp[0] = sin(theta)*sin(z + M_PI/2.0);
1289  flux(0,0) = -((cos(theta)*cos(z)*sin(theta))/r);
1290  flux(1,0) = (cos(theta)*cos(theta)*cos(z))/r;
1291  flux(2,0) = -sin(theta)*sin(z);
1292  flux(0,0) *= -1.;
1293  flux(1,0) *= -1.;
1294  flux(2,0) *= -1.;
1295  }
1296  else if(probAtSphere)
1297  {
1298  //+======================================================================================
1299  // Complete Sphere
1300  REAL r = sqrt(x*x+y*y+z*z);
1301  REAL theta = acos(z/r);
1302  REAL phi = atan2(y,x);
1303 
1304  REAL costheta = cos(theta);
1305  REAL sintheta = sin(theta);
1306 // REAL sin2theta = sin(2.0*theta);
1307 
1308 
1309 
1310  REAL sinphi = sin(phi);
1311  REAL cosphi = cos(phi);
1312  REAL sin2phi = sin(2.0*phi);
1313 
1314  REAL oneminuscosphicosphi = (1.0-cosphi*cosphi);
1315 
1316  REAL npowerofsintheta = 1.0;
1317 
1318  for (int i = 1; i < norder ; i++)
1319  {
1320  npowerofsintheta *= sintheta;
1321  }
1322 
1323  solp[0] = npowerofsintheta*sintheta*oneminuscosphicosphi;
1324 
1325  // Gradient computations
1326  REAL Thetaunitx = cosphi*costheta;
1327  REAL Thetaunity = costheta*sinphi;
1328  REAL Thetaunitz = -sintheta;
1329 
1330  REAL Phiunitx = -sinphi;
1331  REAL Phiunity = cosphi;
1332  REAL Phiunitz = 0.0;
1333 
1334  REAL dfdTheta = (REAL(norder)/r)*costheta*npowerofsintheta*sinphi*sinphi;
1335  REAL dfdPhi = (1.0/r)*npowerofsintheta*sin2phi;
1336 
1337  flux(0,0) = -1.0*(dfdTheta * Thetaunitx + dfdPhi * Phiunitx);
1338  flux(1,0) = -1.0*(dfdTheta * Thetaunity + dfdPhi * Phiunity);
1339  flux(2,0) = -1.0*(dfdTheta * Thetaunitz + dfdPhi * Phiunitz);
1340  }
1341  else
1342  {
1343  std::cout << " Something is wrong " << std::endl;
1344  // One of the previous choices should be true
1345  DebugStop();
1346  }
1347 
1348 }
1349 
1350 void hdivCurvedJCompAppMath::Forcing(const TPZVec<REAL> &pt, TPZVec<STATE> &ff){
1351 
1352  REAL x = pt[0];
1353  REAL y = pt[1];
1354  REAL z = pt[2];
1355  ff.resize(1);
1356 
1357  if(probAtCircle)
1358  {
1359  //+======================================================================================
1360  // Circle
1361  ff[0] = 0.0;
1362  }
1363  else if(probAtCylinder)
1364  {
1365  //+======================================================================================
1366  // Cylinder
1367  REAL r = sqrt(x*x+y*y);
1368  REAL theta = atan2(y,x);
1369  ff[0] = ((1.0 + r*r)*cos(z)*sin(theta))/(r*r);
1370  }
1371  else if(probAtSphere)
1372  {
1373  //+======================================================================================
1374  // Complete Sphere
1375  REAL r = sqrt(x*x+y*y+z*z);
1376  REAL theta = acos(z/r);
1377  REAL phi = atan2(y,x);
1378 
1379  REAL sintheta = sin(theta);
1380  REAL sinphi = sin(phi);
1381  REAL cosphi = cos(phi);
1382  REAL costheta = cos(theta);
1383 
1384  REAL npowerofsintheta = 1.0;
1385  for (int i = 1; i < norder - 1; i++)
1386  {
1387  npowerofsintheta *= sintheta;
1388  }
1389 
1390  ff[0] = -(npowerofsintheta/r*r)*(- REAL(norder) * sintheta*sintheta + cosphi*cosphi*(2.0 + REAL(norder)*sintheta*sintheta)
1391  + (-2.0 + REAL(norder)*REAL(norder)*costheta*costheta)*sinphi*sinphi);
1392  }
1393  else
1394  {
1395  std::cout << " Something is wrong " << std::endl;
1396  // One of the previous choices should be true
1397  DebugStop();
1398  }
1399 
1400 }
1401 
1402 void hdivCurvedJCompAppMath::SolExataH1(const TPZVec<REAL> &pt, TPZVec<STATE> &solp, TPZFMatrix<STATE> &flux){
1403 
1404  solp.resize(1);
1405  solp[0]=0.;
1406  flux.Resize(3, 1);
1407 
1408  // not implemented
1409  DebugStop();
1410 
1411  if(probAtCircle)
1412  {
1413 
1414  }
1415  else if(probAtCylinder)
1416  {
1417 
1418  }
1419  else if(probAtSphere)
1420  {
1421 
1422  }
1423  else
1424  {
1425  std::cout << " Something is wrong " << std::endl;
1426  // One of the previous choices should be true
1427  DebugStop();;
1428  }
1429 
1430 }
1431 
1432 void hdivCurvedJCompAppMath::ForcingH1(const TPZVec<REAL> &pt, TPZVec<STATE> &ff, TPZFMatrix<STATE> &flux)
1433 {
1434  // not implemented
1435  DebugStop();
1436 
1437 
1438  if(probAtCircle)
1439  {
1440 
1441  }
1442  else if(probAtCylinder)
1443  {
1444 
1445  }
1446  else if(probAtSphere)
1447  {
1448 
1449  }
1450  else
1451  {
1452  std::cout << " Something is wrong " << std::endl;
1453  // One of the previous choices should be true
1454  DebugStop();
1455  }
1456 
1457 }
1458 
1459 
1460 void hdivCurvedJCompAppMath::ForcingBC1D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1461  solp.resize(1);
1462  solp[0]=0.;
1463 
1464  REAL x = pt[0];
1465  REAL y = pt[1];
1466  REAL z = pt[2];
1467 
1468  if(probAtCircle)
1469  {
1470  //+======================================================================================
1471  // Circle
1472 
1473  REAL r = sqrt(x*x+y*y+z*z);
1474  REAL theta = atan2(y,x);
1475 
1476  REAL cos2theta = cos(2.0*theta);
1477  REAL sin2theta = sin(2.0*theta);
1478  REAL r4 = r*r*r*r;
1479 
1480  solp[0] = r4*cos2theta*sin2theta;
1481  }
1482  else if(probAtCylinder)
1483  {
1484  //+======================================================================================
1485  // Cylinder
1486 
1487  REAL theta = atan2(y,x);
1488 
1489  solp[0] = sin(theta)*sin(z + M_PI/2.0);
1490  }
1491  else if(probAtSphere)
1492  {
1493  //+======================================================================================
1494  // Complete Sphere
1495  // sobre a casca da esfera -- conferir o raio aqui usado com o da malha geometrica
1496 
1497  REAL theta = atan2(sqrt(x*x+y*y),z);
1498  REAL a = M_PI;
1499 
1500  solp[0] = (a-theta)*sin(theta)*sin(theta);
1501  }
1502  else
1503  {
1504  std::cout << " Something is wrong " << std::endl;
1505  // One of the previous choices should be true
1506  DebugStop();
1507  }
1508 
1509 }
1510 
1511 void hdivCurvedJCompAppMath::ForcingBC2D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1512 
1513  solp.resize(1);
1514  solp[0]=0.;
1515 
1516  REAL x = pt[0];
1517  REAL y = pt[1];
1518  REAL z = pt[2];
1519 
1520  if(probAtCircle)
1521  {
1522  //+======================================================================================
1523  // Circle
1524 
1525  REAL r = sqrt(x*x+y*y+z*z);
1526  REAL theta = atan2(y,x);
1527 
1528  REAL cos2theta = cos(2.0*theta);
1529  REAL sin2theta = sin(2.0*theta);
1530  REAL r4 = r*r*r*r;
1531 
1532  solp[0] = r4*cos2theta*sin2theta;
1533  }
1534  else if(probAtCylinder)
1535  {
1536  //+======================================================================================
1537  // Cylinder
1538 
1539  REAL theta = atan2(y,x);
1540 
1541  solp[0] = sin(theta)*sin(z + M_PI/2.0);
1542  }
1543  else if(probAtSphere)
1544  {
1545  //+======================================================================================
1546  // Complete Sphere
1547  // sobre a casca da esfera -- conferir o raio aqui usado com o da malha geometrica
1548 
1549  REAL theta = atan2(sqrt(x*x+y*y),z);
1550  REAL a = M_PI;
1551 
1552  solp[0] = (a-theta)*sin(theta)*sin(theta);
1553  }
1554  else
1555  {
1556  std::cout << " Something is wrong " << std::endl;
1557  // One of the previous choices should be true
1558  DebugStop();
1559  }
1560 
1561 }
1562 
1563 
1564 void hdivCurvedJCompAppMath::ForcingBC3D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1565 
1566  solp.resize(1);
1567  solp[0]=0.;
1568 
1569  REAL x = pt[0];
1570  REAL y = pt[1];
1571  REAL z = pt[2];
1572 
1573  if(probAtCircle)
1574  {
1575  //+======================================================================================
1576  // Circle
1577 
1578  REAL r = sqrt(x*x+y*y+z*z);
1579  REAL theta = atan2(y,x);
1580 
1581  REAL cos2theta = cos(2.0*theta);
1582  REAL sin2theta = sin(2.0*theta);
1583  REAL r4 = r*r*r*r;
1584 
1585  solp[0] = r4*cos2theta*sin2theta;
1586  }
1587  else if(probAtCylinder)
1588  {
1589  //+======================================================================================
1590  // Cylinder
1591 
1592  REAL theta = atan2(y,x);
1593 
1594  solp[0] = sin(theta)*sin(z + M_PI/2.0);
1595  }
1596  else if(probAtSphere)
1597  {
1598  //+======================================================================================
1599  // Complete Sphere
1600  // sobre a casca da esfera -- conferir o raio aqui usado com o da malha geometrica
1601 
1602  REAL theta = atan2(sqrt(x*x+y*y),z);
1603  REAL a = M_PI;
1604 
1605  solp[0] = (a-theta)*sin(theta)*sin(theta);
1606  }
1607  else
1608  {
1609  std::cout << " Something is wrong " << std::endl;
1610  // One of the previous choices should be true
1611  DebugStop();
1612  }
1613 
1614 }
1615 
1616 
1617 void hdivCurvedJCompAppMath::ForcingBC4D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1618 
1619  solp.resize(1);
1620  solp[0]=0.;
1621 
1622  REAL x = pt[0];
1623  REAL y = pt[1];
1624  REAL z = pt[2];
1625 
1626  if(probAtCircle)
1627  {
1628  //+======================================================================================
1629  // Circle
1630 
1631  REAL r = sqrt(x*x+y*y+z*z);
1632  REAL theta = atan2(y,x);
1633 
1634  REAL cos2theta = cos(2.0*theta);
1635  REAL sin2theta = sin(2.0*theta);
1636  REAL r4 = r*r*r*r;
1637 
1638  solp[0] = r4*cos2theta*sin2theta;
1639  }
1640  else if(probAtCylinder)
1641  {
1642  //+======================================================================================
1643  // Cylinder
1644  REAL theta = atan2(y,x);
1645 
1646  solp[0] = sin(theta)*sin(z + M_PI/2.0);
1647  }
1648  else if(probAtSphere)
1649  {
1650  //+======================================================================================
1651  // Complete Sphere
1652 
1653  REAL theta = atan2(sqrt(x*x+y*y),z);
1654  REAL a = M_PI;
1655 
1656  solp[0] = (a-theta)*sin(theta)*sin(theta);
1657  }
1658  else
1659  {
1660  std::cout << " Something is wrong " << std::endl;
1661  // One of the previous choices should be true
1662  DebugStop();
1663  }
1664 
1665 }
1666 void hdivCurvedJCompAppMath::ForcingBC5D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1667 
1668  REAL x = pt[0];
1669  REAL y = pt[1];
1670  REAL z = pt[2];
1671 
1672 
1673  REAL r = sqrt(x*x+y*y+z*z);
1674  REAL theta = acos(z/r);
1675  REAL phi = atan2(y,x);
1676 
1677 // REAL costheta = cos(theta);
1678  REAL sintheta = sin(theta);
1679  // REAL sin2theta = sin(2.0*theta);
1680 
1681 // REAL sinphi = sin(phi);
1682  REAL cosphi = cos(phi);
1683 // REAL sin2phi = sin(2.0*phi);
1684 
1685  REAL oneminuscosphicosphi = (1.0-cosphi*cosphi);
1686 
1687  REAL npowerofsintheta = 1.0;
1688 
1689  for (int i = 1; i < norder ; i++)
1690  {
1691  npowerofsintheta *= sintheta;
1692  }
1693 
1694  solp[0] = npowerofsintheta*sintheta*oneminuscosphicosphi;
1695 
1696 }
1697 
1698 
1699 
1700 TPZCompMesh *hdivCurvedJCompAppMath::CMeshH1(TPZGeoMesh *gmesh, int pOrder, int dim)
1701 {
1703  TPZMatPoisson3d *material = new TPZMatPoisson3d(fmatId,fDim);
1704  TPZMaterial * mat(material);
1705  material->NStateVariables();
1706 
1707  // //solucao exata
1708  TPZAutoPointer<TPZFunction<STATE> > solexata;
1709  solexata = new TPZDummyFunction<STATE>(SolExataH1, 5);
1710  material->SetForcingFunctionExact(solexata);
1711 
1712  //funcao do lado direito da equacao do problema
1713  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(ForcingH1, 5);
1714  TPZAutoPointer<TPZFunction<STATE> > forcef;
1715  dum->SetPolynomialOrder(20);
1716  forcef = dum;
1717  material->SetForcingFunction(forcef);
1718 
1719  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1720  cmesh->SetDimModel(fDim);
1721  cmesh->SetDefaultOrder(pOrder);
1722 
1723  cmesh->InsertMaterialObject(mat);
1724 
1726  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1727  val2(0,0) = 0.0;
1728  val2(1,0) = 0.0;
1729 
1730  TPZMaterial * BCond1;
1731  TPZMaterial * BCond2;
1732  TPZMaterial * BCond3;
1733  TPZMaterial * BCond4;
1734 
1735 
1736  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1737  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1738  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1739  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1740 
1741  cmesh->InsertMaterialObject(BCond1);
1742  cmesh->InsertMaterialObject(BCond2);
1743  cmesh->InsertMaterialObject(BCond3);
1744  cmesh->InsertMaterialObject(BCond4);
1745 
1746 
1747  cmesh->SetAllCreateFunctionsContinuous();
1748 
1749  //Ajuste da estrutura de dados computacional
1750  cmesh->AutoBuild();
1751 
1752  return cmesh;
1753 
1754 }
1755 
1756 
1757 TPZCompMesh *hdivCurvedJCompAppMath::CMeshFlux(TPZGeoMesh *gmesh, int pOrder, int dim)
1758 {
1759  if(probAtCircle)
1760  {
1762  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
1763 
1764  TPZMaterial * mat(material);
1765  material->NStateVariables();
1766 
1767  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1768 
1770  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1771 
1772  TPZMaterial * BCond1;
1773  TPZMaterial * BCond2;
1774  TPZMaterial * BCond3;
1775  TPZMaterial * BCond4;
1776  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1777  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1778  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1779  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1780 
1781  cmesh->InsertMaterialObject(mat);
1782 
1783  cmesh->SetDimModel(fDim);
1784 
1785  cmesh->SetAllCreateFunctionsHDiv();
1786 
1787 
1788  cmesh->InsertMaterialObject(BCond1);
1789  cmesh->InsertMaterialObject(BCond2);
1790  cmesh->InsertMaterialObject(BCond3);
1791  cmesh->InsertMaterialObject(BCond4);
1792 
1793  cmesh->SetDefaultOrder(pOrder);
1794 
1795  TPZLagrangeMultiplier *matskelet = new TPZLagrangeMultiplier(fmatskeleton, fDim-1, 1);
1796  TPZMaterial * mat2(matskelet);
1797  cmesh->InsertMaterialObject(mat2);
1798 
1799  //Ajuste da estrutura de dados computacional
1800  cmesh->AutoBuild();
1801 
1802  return cmesh;
1803  }
1804  else if(probAtCylinder)
1805  {
1807  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
1808  TPZMaterial * mat(material);
1809  material->NStateVariables();
1810 
1811  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1812 
1814  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1815  TPZMaterial * BCond0;
1816  TPZMaterial * BCond1;
1817  TPZMaterial * BCond2;
1818  TPZMaterial * BCond3;
1819  TPZMaterial * BCond4;
1820  TPZMaterial * BCond5;
1821 
1822  if( dim == 3 ) { BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);}
1823  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1824  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1825  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1826  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1827  if( dim == 3 ) { BCond5 = material->CreateBC(mat, fbc5,fdirichlet, val1, val2);}
1828 
1829  cmesh->InsertMaterialObject(mat);
1830 
1831  cmesh->SetDimModel(dim);
1832 
1833  cmesh->SetAllCreateFunctionsHDiv();
1834 
1835  if( dim == 3 ) { cmesh->InsertMaterialObject(BCond0); }
1836  cmesh->InsertMaterialObject(BCond1);
1837  cmesh->InsertMaterialObject(BCond2);
1838  cmesh->InsertMaterialObject(BCond3);
1839  cmesh->InsertMaterialObject(BCond4);
1840  if( dim == 3 ) { cmesh->InsertMaterialObject(BCond5); }
1841 
1842  cmesh->SetDefaultOrder(pOrder);
1843 
1844  //Ajuste da estrutura de dados computacional
1845  cmesh->AutoBuild();
1846 
1847  return cmesh;
1848  }
1849  else if(probAtSphere)
1850  {
1852  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
1853 
1854  TPZMaterial * mat(material);
1855  material->NStateVariables();
1856 
1857  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1858 
1860  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1861  TPZMaterial * BCond0;
1862  TPZMaterial * BCond1;
1863 
1864  cmesh->InsertMaterialObject(mat);
1865  cmesh->SetDimModel(dim);
1866  cmesh->SetAllCreateFunctionsHDiv();
1867 
1868  BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);
1869  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1870  cmesh->InsertMaterialObject(BCond0);
1871  cmesh->InsertMaterialObject(BCond1);
1872 
1873  cmesh->SetDefaultOrder(pOrder);
1874 
1875  TPZLagrangeMultiplier *matskelet = new TPZLagrangeMultiplier(fmatskeleton, dim-1, 1);
1876  TPZMaterial * mat2(matskelet);
1877  cmesh->InsertMaterialObject(mat2);
1878 
1879  //Ajuste da estrutura de dados computacional
1880  cmesh->AutoBuild();
1881 
1882  this->SetupDisconnectedHdivboud(fbc0,fbc1,cmesh);
1883 
1884  return cmesh;
1885  }
1886  else
1887  {
1888  std::cout << " Something is wrong " << std::endl;
1889  // One of the previous choices should be true
1890  DebugStop();
1891  }
1892  return new TPZCompMesh(gmesh);
1893 }
1894 
1895 void hdivCurvedJCompAppMath::SetupDisconnectedHdivboud(const int left,const int rigth, TPZCompMesh * cmesh)
1896 {
1897 
1898  int64_t ncels = cmesh->NElements();
1899  for (int64_t icel = 0; icel < ncels; icel++) {
1900  TPZCompEl * cel = cmesh->Element(icel);
1901  if (!cel) {
1902  DebugStop();
1903  }
1904 
1905  if (cel->Dimension() != 1) {
1906  continue;
1907  }
1908 
1909  if (cel->Reference()->MaterialId() == fbc0)
1910  {
1911  TPZConnect & connect = cel->Connect(0);
1912  int64_t newindex = cmesh->AllocateNewConnect(connect);
1913  TPZGeoEl * gel = cel->Reference();
1914  if (!gel) {
1915  DebugStop();
1916  }
1917 
1918  TPZGeoElSide gelside(gel,gel->NSides()-1);
1919 
1920  TPZStack<TPZCompElSide > celstack;
1921  gelside.EqualLevelCompElementList(celstack, 1, 0);
1922 
1923  int64_t sz = celstack.size();
1924  if (sz == 0) {
1925  DebugStop();
1926  }
1927 
1928  for (int64_t iside=0; iside < sz; iside++) {
1929  TPZCompElSide cels = celstack[iside];
1930 
1931  if(cels.Element()->Dimension() == 2)
1932  {
1933  TPZInterpolationSpace * InterpElside = dynamic_cast<TPZInterpolationSpace *>(cels.Element());
1934  int nsideconnects = InterpElside->NSideConnects(cels.Side());
1935  if(nsideconnects != 1 )
1936  {
1937  DebugStop();
1938  }
1939  int localindex = InterpElside->SideConnectLocId(0, cels.Side());
1940  InterpElside->SetConnectIndex(localindex, newindex);
1941 // int orientside = InterpElside->GetSideOrient(cels.Side());
1942  InterpElside->SetSideOrient(cels.Side(),1);
1943 
1944  cel->SetConnectIndex(0, newindex);
1945  TPZInterpolationSpace * celbound = dynamic_cast<TPZInterpolationSpace *>(cel);
1946  celbound->SetSideOrient(2,1);
1947  cmesh->ExpandSolution();
1948 // cmesh->CleanUpUnconnectedNodes();
1949  break;
1950  }
1951  }
1952  }
1953  }
1954 
1955 }
1956 
1957 TPZCompMesh *hdivCurvedJCompAppMath::CMeshPressure(TPZGeoMesh *gmesh, int pOrder, int dim)
1958 {
1959 
1960  if(probAtCircle)
1961  {
1962 
1964  TPZMatPoisson3d *material = new TPZMatPoisson3d(fmatId,fDim);
1965  material->NStateVariables();
1966 
1967  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1968  cmesh->SetDimModel(fDim);
1969  TPZMaterial * mat(material);
1970  cmesh->InsertMaterialObject(mat);
1971 
1973  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1974 
1975  cmesh->SetDefaultOrder(pOrder);
1976 
1977  bool h1function = true;// com esqueleto precisa disso
1978  if(h1function){
1979  cmesh->SetAllCreateFunctionsContinuous();
1980  cmesh->ApproxSpace().CreateDisconnectedElements(true);
1981  }
1982  else{
1983  cmesh->SetAllCreateFunctionsDiscontinuous();
1984  }
1985 
1986  //Ajuste da estrutura de dados computacional
1987  cmesh->AutoBuild();
1988 
1989 
1990  int ncon = cmesh->NConnects();
1991  for(int i=0; i<ncon; i++)
1992  {
1993  TPZConnect &newnod = cmesh->ConnectVec()[i];
1994  //newnod.SetPressure(true);
1995  newnod.SetLagrangeMultiplier(1);
1996  }
1997 
1998  if(!h1function)
1999  {
2000 
2001  int nel = cmesh->NElements();
2002  for(int i=0; i<nel; i++){
2003  TPZCompEl *cel = cmesh->ElementVec()[i];
2004  TPZCompElDisc *celdisc = dynamic_cast<TPZCompElDisc *>(cel);
2005  celdisc->SetConstC(1.);
2006  celdisc->SetCenterPoint(0, 0.);
2007  celdisc->SetCenterPoint(1, 0.);
2008  celdisc->SetCenterPoint(2, 0.);
2009  celdisc->SetTrueUseQsiEta();
2010  //celdisc->SetFalseUseQsiEta();
2011 
2012  if(celdisc && celdisc->Reference()->Dimension() == cmesh->Dimension())
2013  {
2014  if(ftriang==true) celdisc->SetTotalOrderShape();
2015  else celdisc->SetTensorialShape();
2016  }
2017 
2018  }
2019  }
2020 
2021 #ifdef PZDEBUG
2022  int ncel = cmesh->NElements();
2023  for(int i =0; i<ncel; i++){
2024  TPZCompEl * compEl = cmesh->ElementVec()[i];
2025  if(!compEl) continue;
2026  TPZInterfaceElement * facel = dynamic_cast<TPZInterfaceElement *>(compEl);
2027  if(facel)DebugStop();
2028 
2029  }
2030 #endif
2031 
2032 // int nel = cmesh->NElements();
2033 // for(int i=0; i<nel; i++){
2034 // TPZCompEl *cel = cmesh->ElementVec()[i];
2035 //
2036 //
2037 // }
2038  return cmesh;
2039  }
2040  else if(probAtCylinder)
2041  {
2043  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
2044  material->NStateVariables();
2045 
2046  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
2047  cmesh->SetDimModel(dim);
2048  TPZMaterial * mat(material);
2049  cmesh->InsertMaterialObject(mat);
2050 
2052  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
2053 
2054  cmesh->SetDefaultOrder(pOrder);
2055 
2056  bool h1function = true;//false em triangulo
2057  if(h1function){
2058  cmesh->SetAllCreateFunctionsContinuous();
2059  cmesh->ApproxSpace().CreateDisconnectedElements(true);
2060  }
2061  else{
2062  cmesh->SetAllCreateFunctionsDiscontinuous();
2063  }
2064 
2065 
2066  //Ajuste da estrutura de dados computacional
2067  cmesh->AutoBuild();
2068 
2069 
2070  int ncon = cmesh->NConnects();
2071  for(int i=0; i<ncon; i++)
2072  {
2073  TPZConnect &newnod = cmesh->ConnectVec()[i];
2074  //newnod.SetPressure(true);
2075  newnod.SetLagrangeMultiplier(1);
2076  }
2077 
2078  if(!h1function)
2079  {
2080 
2081  int nel = cmesh->NElements();
2082  for(int i=0; i<nel; i++){
2083  TPZCompEl *cel = cmesh->ElementVec()[i];
2084  TPZCompElDisc *celdisc = dynamic_cast<TPZCompElDisc *>(cel);
2085  celdisc->SetConstC(1.);
2086  celdisc->SetCenterPoint(0, 0.);
2087  celdisc->SetCenterPoint(1, 0.);
2088  celdisc->SetCenterPoint(2, 0.);
2089  celdisc->SetTrueUseQsiEta();
2090  //celdisc->SetFalseUseQsiEta();
2091 
2092  if(celdisc && celdisc->Reference()->Dimension() == cmesh->Dimension())
2093  {
2094  if(ftriang) celdisc->SetTotalOrderShape();
2095  else celdisc->SetTensorialShape();
2096  }
2097 
2098  }
2099  }
2100 
2101 #ifdef PZDEBUG
2102  int ncel = cmesh->NElements();
2103  for(int i =0; i<ncel; i++){
2104  TPZCompEl * compEl = cmesh->ElementVec()[i];
2105  if(!compEl) continue;
2106  TPZInterfaceElement * facel = dynamic_cast<TPZInterfaceElement *>(compEl);
2107  if(facel)DebugStop();
2108 
2109  }
2110 #endif
2111 
2112  return cmesh;
2113  }
2114  else if(probAtSphere)
2115  {
2116  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
2117 
2118  material->NStateVariables();
2119 
2120  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
2121  cmesh->SetDimModel(dim);
2122  TPZMaterial * mat(material);
2123  cmesh->InsertMaterialObject(mat);
2124 
2126  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
2127 
2128  cmesh->SetDefaultOrder(pOrder);
2129  //cmesh->SetDimModel(dim);
2130 
2131  bool h1function = true;//(espaco sera sempre assim quando estiver usando elementos wrap)
2132  if(h1function){
2133  cmesh->SetAllCreateFunctionsContinuous();
2134  cmesh->ApproxSpace().CreateDisconnectedElements(true);
2135  }
2136  else{
2137  cmesh->SetAllCreateFunctionsDiscontinuous();
2138  }
2139 
2140 
2141  //Ajuste da estrutura de dados computacional
2142  cmesh->AutoBuild();
2143 
2144 
2145  int ncon = cmesh->NConnects();
2146  for(int i=0; i<ncon; i++)
2147  {
2148  TPZConnect &newnod = cmesh->ConnectVec()[i];
2149  //newnod.SetPressure(true);
2150  newnod.SetLagrangeMultiplier(1);
2151  }
2152 
2153  if(!h1function)
2154  {
2155 
2156  int nel = cmesh->NElements();
2157  for(int i=0; i<nel; i++){
2158  TPZCompEl *cel = cmesh->ElementVec()[i];
2159  TPZCompElDisc *celdisc = dynamic_cast<TPZCompElDisc *>(cel);
2160  celdisc->SetConstC(1.);
2161  celdisc->SetCenterPoint(0, 0.);
2162  celdisc->SetCenterPoint(1, 0.);
2163  celdisc->SetCenterPoint(2, 0.);
2164  celdisc->SetTrueUseQsiEta();
2165  //celdisc->SetFalseUseQsiEta();
2166 
2167 
2168  if(celdisc && celdisc->Reference()->Dimension() == cmesh->Dimension())
2169  {
2170  if(ftriang==true) celdisc->SetTotalOrderShape();
2171  else celdisc->SetTensorialShape();
2172  }
2173 
2174  }
2175  }
2176 
2177  //#endif
2178  return cmesh;
2179  }
2180  else
2181  {
2182  std::cout << " Something is wrong " << std::endl;
2183  // One of the previous choices should be true
2184  DebugStop();
2185  }
2186  return new TPZCompMesh;
2187 }
2188 
2189 TPZCompMesh *hdivCurvedJCompAppMath::CMeshMixed(TPZGeoMesh * gmesh, TPZVec<TPZCompMesh *> meshvec)
2190 {
2191 
2192  if(probAtCircle)
2193  {
2194 
2195  //Creating computational mesh for multiphysic elements
2196  gmesh->ResetReference();
2197  TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
2198 
2199  //criando material
2200  int dim = gmesh->Dimension();
2201  bool intface;
2202 // TPZMatPoissonD3 *material = new TPZMatPoissonD3(fmatId,dim); intface = true; // nesse material tem que ser true
2203  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,dim); intface = true; // nesse material tem que ser true
2204  //TPZMixedPoisson *material = new TPZMixedPoisson(matId,dim); intface = false; // nesse material tem que ser false
2205 
2206  //incluindo os dados do problema
2207 
2208  //incluindo os dados do problema
2209  TPZFNMatrix<2,REAL> PermTensor(dim,dim,0.);
2210  TPZFNMatrix<2,REAL> InvPermTensor(dim,dim,0.);
2211 
2212  // tensor de permutacao
2213  TPZFNMatrix<9,REAL> TP(dim,dim,0.0);
2214  TPZFNMatrix<9,REAL> InvTP(dim,dim,0.0);
2215 
2216  // Hard coded
2217  for (int id = 0; id < dim; id++){
2218  TP(id,id) = 1.0;
2219  InvTP(id,id) = 1.0;
2220  }
2221 
2222  PermTensor = TP;
2223  InvPermTensor = InvTP;
2224 
2225  material->SetPermeabilityTensor(PermTensor, InvPermTensor);
2226 
2227  //solucao exata
2228  TPZAutoPointer<TPZFunction<STATE> > solexata;
2229 
2230  solexata = new TPZDummyFunction<STATE>(SolExata, 5);
2231  material->SetForcingFunctionExact(solexata);
2232  mphysics->SetDimModel(dim);
2233  //funcao do lado direito da equacao do problema
2234  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(Forcing, 5);
2235  TPZAutoPointer<TPZFunction<STATE> > forcef;
2236  dum->SetPolynomialOrder(20);
2237  forcef = dum;
2238  material->SetForcingFunction(forcef);
2239 
2240  //inserindo o material na malha computacional
2241  TPZMaterial *mat(material);
2242  mphysics->InsertMaterialObject(mat);
2243 
2244 
2245  //Criando condicoes de contorno
2246  TPZMaterial * BCond0 = NULL;
2247  TPZMaterial * BCond1;
2248  TPZMaterial * BCond2 = NULL;
2249  TPZMaterial * BCond3 = NULL;
2250  TPZMaterial * BCond4 = NULL;
2251  TPZMaterial * BCond5 = NULL;
2252 
2253  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
2254 
2255 
2256  val2(0,0) = 0.0;
2257  val2(1,0) = 0.0;
2258  TPZAutoPointer<TPZFunction<STATE> > FBCond1 = new TPZDummyFunction<STATE>(ForcingBC1D, 5);
2259  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
2260  BCond1->SetForcingFunction(FBCond1);
2261  // BCond1 = material->CreateBC(mat, bc1,neumann, val1, val2);
2262 
2263  mphysics->SetAllCreateFunctionsMultiphysicElem();
2264  if( dim == 3 ) { mphysics->InsertMaterialObject(BCond0); }
2265  mphysics->InsertMaterialObject(BCond1);
2266  mphysics->InsertMaterialObject(BCond2);
2267  mphysics->InsertMaterialObject(BCond3);
2268  mphysics->InsertMaterialObject(BCond4);
2269  if( dim == 3 ) { mphysics->InsertMaterialObject(BCond5); }
2270 
2271  //Fazendo auto build
2272  mphysics->AutoBuild();
2273  mphysics->AdjustBoundaryElements();
2274  mphysics->CleanUpUnconnectedNodes();
2275 
2276  // Creating multiphysic elements into mphysics computational mesh
2277  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2278  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
2279  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
2280 
2281  mphysics->Reference()->ResetReference();
2282  mphysics->LoadReferences();
2283 
2284  // Creation of interface elements
2285  if (intface)
2286  {
2287  int nel = mphysics->ElementVec().NElements();
2288  for(int el = 0; el < nel; el++)
2289  {
2290  TPZCompEl * compEl = mphysics->ElementVec()[el];
2291  if(!compEl) continue;
2292  int index = compEl ->Index();
2293  if(compEl->Dimension() == mphysics->Dimension())
2294  {
2295  TPZMultiphysicsElement * InterpEl = dynamic_cast<TPZMultiphysicsElement *>(mphysics->ElementVec()[index]);
2296  if(!InterpEl) continue;
2297  InterpEl->CreateInterfaces();
2298  }
2299  }
2300 
2301  }
2302 
2303  return mphysics;
2304  }
2305  else if(probAtCylinder)
2306  {
2307  bool condensacaoestatica = true;
2308 
2309  //Creating computational mesh for multiphysic elements
2310  gmesh->ResetReference();
2311  TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
2312 
2313  //criando material
2314  int dim = gmesh->Dimension();
2315 
2316  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,dim);
2317 
2318  //incluindo os dados do problema
2319  TPZFNMatrix<9,REAL> PermTensor(3,3,0.);
2320  TPZFNMatrix<9,REAL> InvPermTensor(3,3,0.);
2321 
2322 
2323  // tensor de permutacao
2324  TPZFNMatrix<9,REAL> TP(3,3,0.0);
2325  TPZFNMatrix<9,REAL> InvTP(3,3,0.0);
2326 
2327  // Hard coded
2328  for (int id = 0; id < 3; id++){
2329  TP(id,id) = 1.0;
2330  InvTP(id,id) = 1.0;
2331  }
2332 
2333  PermTensor = TP;
2334  InvPermTensor = InvTP;
2335 
2336  material->SetPermeabilityTensor(PermTensor, InvPermTensor);
2337 
2338  //solucao exata
2339  TPZAutoPointer<TPZFunction<STATE> > solexata;
2340 
2341  solexata = new TPZDummyFunction<STATE>(SolExata, 5);
2342  material->SetForcingFunctionExact(solexata);
2343  mphysics->SetDimModel(dim);
2344  //funcao do lado direito da equacao do problema
2345  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(Forcing, 5);
2346  TPZAutoPointer<TPZFunction<STATE> > forcef;
2347  dum->SetPolynomialOrder(1);
2348  forcef = dum;
2349  material->SetForcingFunction(forcef);
2350 
2351  //inserindo o material na malha computacional
2352  TPZMaterial *mat(material);
2353  mphysics->InsertMaterialObject(mat);
2354 
2355 
2356  //Criando condicoes de contorno
2357  TPZMaterial * BCond0 = NULL;
2358  TPZMaterial * BCond1;
2359  TPZMaterial * BCond2;
2360  TPZMaterial * BCond3;
2361  TPZMaterial * BCond4;
2362  TPZMaterial * BCond5 = NULL;
2363 
2364  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
2365 
2366  val2(0,0) = 0.0;
2367  val2(1,0) = 0.0;
2368  TPZAutoPointer<TPZFunction<STATE> > FBCond1 = new TPZDummyFunction<STATE>(ForcingBC1D, 5);
2369  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
2370  BCond1->SetForcingFunction(FBCond1);
2371  // BCond1 = material->CreateBC(mat, bc1,neumann, val1, val2);
2372 
2373  val2(0,0) = 0.0;
2374  val2(1,0) = 0.0;
2375  TPZAutoPointer<TPZFunction<STATE> > FBCond2 = new TPZDummyFunction<STATE>(ForcingBC2D, 5);
2376  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
2377  BCond2->SetForcingFunction(FBCond2);
2378  // TPZAutoPointer<TPZFunction<STATE> > FBCond2 = new TPZDummyFunction<STATE>(ForcingBC2N, 5);
2379  // BCond2 = material->CreateBC(mat, bc2,neumann, val1, val2);
2380  // BCond2->SetForcingFunction(FBCond2);
2381 
2382  val2(0,0) = 0.0;
2383  val2(1,0) = 0.0;
2384  TPZAutoPointer<TPZFunction<STATE> > FBCond3 = new TPZDummyFunction<STATE>(ForcingBC3D, 5);
2385  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
2386  BCond3->SetForcingFunction(FBCond3);
2387  // BCond3 = material->CreateBC(mat, bc3,neumann, val1, val2);
2388 
2389  val2(0,0) = 0.0;
2390  val2(1,0) = 0.0;
2391  TPZAutoPointer<TPZFunction<STATE> > FBCond4 = new TPZDummyFunction<STATE>(ForcingBC4D, 5);
2392  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
2393  BCond4->SetForcingFunction(FBCond4);
2394 
2395 
2396  mphysics->SetAllCreateFunctionsMultiphysicElem();
2397  if( dim == 3 ) { mphysics->InsertMaterialObject(BCond0); }
2398  mphysics->InsertMaterialObject(BCond1);
2399  mphysics->InsertMaterialObject(BCond2);
2400  mphysics->InsertMaterialObject(BCond3);
2401  mphysics->InsertMaterialObject(BCond4);
2402  if( dim == 3 ) { mphysics->InsertMaterialObject(BCond5); }
2403 
2404  //Fazendo auto build
2405  mphysics->AutoBuild();
2406  mphysics->AdjustBoundaryElements();
2407  mphysics->CleanUpUnconnectedNodes();
2408 
2409 
2410  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2411 
2412  if (condensacaoestatica)
2413  {
2414  //Creating multiphysic elements containing skeletal elements.
2415  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2416  mphysics->Reference()->ResetReference();
2417  mphysics->LoadReferences();
2418 
2419  int64_t nel = mphysics->ElementVec().NElements();
2420 
2421  std::map<int64_t, int64_t> bctoel, eltowrap;
2422  for (int64_t el=0; el<nel; el++) {
2423  TPZCompEl *cel = mphysics->Element(el);
2424  TPZGeoEl *gel = cel->Reference();
2425  int matid = gel->MaterialId();
2426  if (matid < 0) {
2427  TPZGeoElSide gelside(gel,gel->NSides()-1);
2428  TPZGeoElSide neighbour = gelside.Neighbour();
2429  while (neighbour != gelside) {
2430  if (neighbour.Element()->Dimension() == dim && neighbour.Element()->Reference()) {
2431  // got you!!
2432  bctoel[el] = neighbour.Element()->Reference()->Index();
2433  break;
2434  }
2435  neighbour = neighbour.Neighbour();
2436  }
2437  if (neighbour == gelside) {
2438  DebugStop();
2439  }
2440  }
2441  }
2442 
2443  TPZStack< TPZStack< TPZMultiphysicsElement *,7> > wrapEl;
2444  for(int64_t el = 0; el < nel; el++)
2445  {
2446  TPZMultiphysicsElement *mfcel = dynamic_cast<TPZMultiphysicsElement *>(mphysics->Element(el));
2447  if(mfcel->Dimension()==dim) TPZBuildMultiphysicsMesh::AddWrap(mfcel, fmatId, wrapEl);//criei elementos com o mesmo matId interno, portanto nao preciso criar elemento de contorno ou outro material do tipo TPZLagrangeMultiplier
2448  }
2449 
2450  for (int64_t el =0; el < wrapEl.size(); el++) {
2451  TPZCompEl *cel = wrapEl[el][0];
2452  int64_t index = cel->Index();
2453  eltowrap[index] = el;
2454  }
2455 
2456  meshvec[0]->CleanUpUnconnectedNodes();
2457  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
2458  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
2459 
2460  std::map<int64_t, int64_t>::iterator it;
2461  for (it = bctoel.begin(); it != bctoel.end(); it++) {
2462  int64_t bcindex = it->first;
2463  int64_t elindex = it->second;
2464  if (eltowrap.find(elindex) == eltowrap.end()) {
2465  DebugStop();
2466  }
2467  int64_t wrapindex = eltowrap[elindex];
2468  TPZCompEl *bcel = mphysics->Element(bcindex);
2469  TPZMultiphysicsElement *bcmf = dynamic_cast<TPZMultiphysicsElement *>(bcel);
2470  if (!bcmf) {
2471  DebugStop();
2472  }
2473  wrapEl[wrapindex].Push(bcmf);
2474 
2475  }
2476 
2477  //------- Create and add group elements -------
2478  int64_t index, nenvel;
2479  nenvel = wrapEl.NElements();
2480  TPZStack<TPZElementGroup *> elgroups;
2481  for(int64_t ienv=0; ienv<nenvel; ienv++){
2482  TPZElementGroup *elgr = new TPZElementGroup(*wrapEl[ienv][0]->Mesh(),index);
2483  elgroups.Push(elgr);
2484  nel = wrapEl[ienv].NElements();
2485  for(int jel=0; jel<nel; jel++){
2486  elgr->AddElement(wrapEl[ienv][jel]);
2487  }
2488  }
2489 
2490  mphysics->ComputeNodElCon();
2491  // create condensed elements
2492  // increase the NumElConnected of one pressure connects in order to prevent condensation
2493  for (int64_t ienv=0; ienv<nenvel; ienv++) {
2494  TPZElementGroup *elgr = elgroups[ienv];
2495  int nc = elgr->NConnects();
2496  for (int ic=0; ic<nc; ic++) {
2497  TPZConnect &c = elgr->Connect(ic);
2498  if (c.LagrangeMultiplier() > 0) {
2499  c.IncrementElConnected();
2500  break;
2501  }
2502  }
2503 
2504  new TPZCondensedCompEl(elgr);
2505  }
2506 
2507  mphysics->CleanUpUnconnectedNodes();
2508  mphysics->ExpandSolution();
2509  }
2510  else
2511  {
2512  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2513  mphysics->Reference()->ResetReference();
2514  mphysics->LoadReferences();
2515 
2516  // TPZMaterial * skeletonEl = material->CreateBC(mat, matskeleton, 3, val1, val2);
2517  // mphysics->InsertMaterialObject(skeletonEl);
2518 
2519  // TPZLagrangeMultiplier *matskelet = new TPZLagrangeMultiplier(matskeleton, dim-1, 1);
2520  // TPZMaterial * mat2(matskelet);
2521  // mphysics->InsertMaterialObject(mat2);
2522 
2523  int nel = mphysics->ElementVec().NElements();
2524  TPZStack< TPZStack< TPZMultiphysicsElement *,7> > wrapEl;
2525  for(int el = 0; el < nel; el++)
2526  {
2527  TPZMultiphysicsElement *mfcel = dynamic_cast<TPZMultiphysicsElement *>(mphysics->Element(el));
2528  if(mfcel->Dimension()==dim) TPZBuildMultiphysicsMesh::AddWrap(mfcel, fmatId, wrapEl);//criei elementos com o mesmo matId interno, portanto nao preciso criar elemento de contorno ou outro material do tipo TPZLagrangeMultiplier
2529  }
2530 
2531  meshvec[0]->CleanUpUnconnectedNodes();
2532  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
2533  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
2534 
2535  // //------- Create and add group elements -------
2536  // int64_t index, nenvel;
2537  // nenvel = wrapEl.NElements();
2538  // for(int ienv=0; ienv<nenvel; ienv++){
2539  // TPZElementGroup *elgr = new TPZElementGroup(*wrapEl[ienv][0]->Mesh(),index);
2540  // nel = wrapEl[ienv].NElements();
2541  // for(int jel=0; jel<nel; jel++){
2542  // elgr->AddElement(wrapEl[ienv][jel]);
2543  // }
2544  // }
2545 
2546  }
2547 
2548 
2549  return mphysics;
2550  }
2551  else if(probAtSphere)
2552  {
2553  bool condensacaoestatica = true; //false; //
2554 
2555  //Creating computational mesh for multiphysic elements
2556  gmesh->ResetReference();
2557  TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
2558 
2559  //criando material
2560  int dim = gmesh->Dimension();
2561 
2562  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,dim);
2563 
2564  //incluindo os dados do problema
2565  TPZFNMatrix<9,REAL> PermTensor(3,3,0.);
2566  TPZFNMatrix<9,REAL> InvPermTensor(3,3,0.);
2567 
2568 
2569  // tensor de permutacao
2570  // Hard coded
2571  for (int id = 0; id < 3; id++){
2572  PermTensor(id,id) = 1.0;
2573  InvPermTensor(id,id) = 1.0;
2574  }
2575 
2576  material->SetPermeabilityTensor(PermTensor, InvPermTensor);
2577 
2578  //solucao exata
2579  TPZAutoPointer<TPZFunction<STATE> > solexata;
2580 
2581  solexata = new TPZDummyFunction<STATE>(SolExata, 5);
2582  material->SetForcingFunctionExact(solexata);
2583  mphysics->SetDimModel(dim);
2584  //funcao do lado direito da equacao do problema
2585  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(Forcing, 5);
2586  TPZAutoPointer<TPZFunction<STATE> > forcef;
2587  dum->SetPolynomialOrder(20);
2588  forcef = dum;
2589  material->SetForcingFunction(forcef);
2590 
2591  //inserindo o material na malha computacional
2592  TPZMaterial *mat(material);
2593  mphysics->InsertMaterialObject(mat);
2594 
2595 
2596  //Criando condicoes de contorno
2597  TPZMaterial * BCond0;
2598  TPZMaterial * BCond1;
2599 
2600  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
2601  val2(0,0) = 0.0;
2602  val2(1,0) = 0.0;
2603  TPZDummyFunction<STATE> * FBCond0D = new TPZDummyFunction<STATE>(ForcingBC5D, 5);
2604  //FBCond0D->SetPolynomialOrder(20);
2605  TPZAutoPointer<TPZFunction<STATE> > FBCond0 = FBCond0D;
2606  BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);
2607  BCond0->SetForcingFunction(FBCond0);
2608 
2609  val2(0,0) = 0.0;
2610  val2(1,0) = 0.0;
2611  TPZDummyFunction<STATE> * FBCond1D = new TPZDummyFunction<STATE>(ForcingBC5D, 5);
2612  //FBCond0D->SetPolynomialOrder(20);
2613  TPZAutoPointer<TPZFunction<STATE> > FBCond1 = FBCond1D;
2614  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
2615  BCond1->SetForcingFunction(FBCond1);
2616 
2617  mphysics->InsertMaterialObject(BCond0);
2618  mphysics->InsertMaterialObject(BCond1);
2619 
2620 
2621  mphysics->SetAllCreateFunctionsMultiphysicElem();
2622 
2623  //Fazendo auto build
2624  mphysics->AutoBuild();
2625  mphysics->AdjustBoundaryElements();
2626  mphysics->CleanUpUnconnectedNodes();
2627 
2628  if (condensacaoestatica)
2629  {
2630  //Creating multiphysic elements containing skeletal elements.
2631  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2632  mphysics->Reference()->ResetReference();
2633  mphysics->LoadReferences();
2634 
2635  int64_t nel = mphysics->ElementVec().NElements();
2636 
2637  std::map<int64_t, int64_t> bctoel, eltowrap;
2638  for (int64_t el=0; el<nel; el++) {
2639  TPZCompEl *cel = mphysics->Element(el);
2640  TPZGeoEl *gel = cel->Reference();
2641  int matid = gel->MaterialId();
2642  if (matid < 0) {
2643  TPZGeoElSide gelside(gel,gel->NSides()-1);
2644  TPZGeoElSide neighbour = gelside.Neighbour();
2645  while (neighbour != gelside) {
2646  if (neighbour.Element()->Dimension() == dim && neighbour.Element()->Reference()) {
2647  // got you!!
2648  bctoel[el] = neighbour.Element()->Reference()->Index();
2649  break;
2650  }
2651  neighbour = neighbour.Neighbour();
2652  }
2653  if (neighbour == gelside) {
2654  DebugStop();
2655  }
2656  }
2657  }
2658 
2659  TPZStack< TPZStack< TPZMultiphysicsElement *,7> > wrapEl;
2660  for(int64_t el = 0; el < nel; el++)
2661  {
2662  TPZMultiphysicsElement *mfcel = dynamic_cast<TPZMultiphysicsElement *>(mphysics->Element(el));
2663  if(mfcel->Dimension()==dim) TPZBuildMultiphysicsMesh::AddWrap(mfcel, fmatId, wrapEl);//criei elementos com o mesmo matId interno, portanto nao preciso criar elemento de contorno ou outro material do tipo TPZLagrangeMultiplier
2664  }
2665 
2666  for (int64_t el =0; el < wrapEl.size(); el++) {
2667  TPZCompEl *cel = wrapEl[el][0];
2668  int64_t index = cel->Index();
2669  eltowrap[index] = el;
2670  }
2671 
2672  meshvec[0]->CleanUpUnconnectedNodes();
2673  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
2674  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
2675 
2676  std::map<int64_t, int64_t>::iterator it;
2677  for (it = bctoel.begin(); it != bctoel.end(); it++) {
2678  int64_t bcindex = it->first;
2679  int64_t elindex = it->second;
2680  if (eltowrap.find(elindex) == eltowrap.end()) {
2681  DebugStop();
2682  }
2683  int64_t wrapindex = eltowrap[elindex];
2684  TPZCompEl *bcel = mphysics->Element(bcindex);
2685  TPZMultiphysicsElement *bcmf = dynamic_cast<TPZMultiphysicsElement *>(bcel);
2686  if (!bcmf) {
2687  DebugStop();
2688  }
2689  wrapEl[wrapindex].Push(bcmf);
2690 
2691  }
2692 
2693  //------- Create and add group elements -------
2694  int64_t index, nenvel;
2695  nenvel = wrapEl.NElements();
2696  TPZStack<TPZElementGroup *> elgroups;
2697  for(int64_t ienv=0; ienv<nenvel; ienv++){
2698  TPZElementGroup *elgr = new TPZElementGroup(*wrapEl[ienv][0]->Mesh(),index);
2699  elgroups.Push(elgr);
2700  nel = wrapEl[ienv].NElements();
2701  for(int jel=0; jel<nel; jel++){
2702  elgr->AddElement(wrapEl[ienv][jel]);
2703  }
2704  }
2705 
2706  mphysics->ComputeNodElCon();
2707  // create condensed elements
2708  // increase the NumElConnected of one pressure connects in order to prevent condensation
2709  for (int64_t ienv=0; ienv<nenvel; ienv++) {
2710  TPZElementGroup *elgr = elgroups[ienv];
2711  int nc = elgr->NConnects();
2712  for (int ic=0; ic<nc; ic++) {
2713  TPZConnect &c = elgr->Connect(ic);
2714  if (c.LagrangeMultiplier() > 0) {
2715  c.IncrementElConnected();
2716  break;
2717  }
2718  }
2719 
2720  new TPZCondensedCompEl(elgr);
2721  }
2722 
2723  mphysics->CleanUpUnconnectedNodes();
2724  mphysics->ExpandSolution();
2725  }
2726  else
2727  {
2728  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
2729  mphysics->Reference()->ResetReference();
2730  mphysics->LoadReferences();
2731 
2732  int nel = mphysics->ElementVec().NElements();
2733  TPZStack< TPZStack< TPZMultiphysicsElement *,7> > wrapEl;
2734  for(int el = 0; el < nel; el++)
2735  {
2736  TPZMultiphysicsElement *mfcel = dynamic_cast<TPZMultiphysicsElement *>(mphysics->Element(el));
2737  if(mfcel->Dimension()==dim) TPZBuildMultiphysicsMesh::AddWrap(mfcel, fmatId, wrapEl);//criei elementos com o mesmo matId interno, portanto nao preciso criar elemento de contorno ou outro material do tipo TPZLagrangeMultiplier
2738  }
2739 
2740  meshvec[0]->CleanUpUnconnectedNodes();
2741  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
2742  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
2743 
2744  }
2745 
2746 
2747 
2748  mphysics->ComputeNodElCon();
2749  mphysics->CleanUpUnconnectedNodes();
2750  mphysics->ExpandSolution();
2751 
2752  return mphysics;
2753  }
2754  else
2755  {
2756  std::cout << " Something is wrong " << std::endl;
2757  // One of the previous choices should be true
2758  DebugStop();
2759  }
2760  return new TPZCompMesh;
2761 }
2762 
2763 
2764 void hdivCurvedJCompAppMath::ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
2765 {
2766 
2767  int64_t nel = l2mesh->NElements();
2768  int dim = l2mesh->Dimension();
2769  TPZManVector<STATE,10> globalerrors(10,0.);
2770  for (int64_t el=0; el<nel; el++) {
2771  TPZCompEl *cel = l2mesh->ElementVec()[el];
2772  if (!cel) {
2773  continue;
2774  }
2775  TPZGeoEl *gel = cel->Reference();
2776  if (!gel || gel->Dimension() != dim) {
2777  continue;
2778  }
2779  TPZManVector<REAL,10> elerror(10,0.);
2780  elerror.Fill(0.);
2781  cel->EvaluateError(SolExataH1, elerror, 0);
2782 
2783  int nerr = elerror.size();
2784  globalerrors.resize(nerr);
2785 
2786  for (int i=0; i<nerr; i++) {
2787  globalerrors[i] += elerror[i]*elerror[i];
2788  }
2789  }
2790  // sequence of storage
2791  // for each
2792  errors.PutVal(pos, 0, p); // polinomial order
2793  errors.PutVal(pos, 1, ndiv); // refinement
2794  errors.PutVal(pos, 2, sqrt(globalerrors[1])); // pressure
2795  errors.PutVal(pos, 3, sqrt(globalerrors[2])); // flux
2796 
2797 }
2798 
2799 void hdivCurvedJCompAppMath::ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, std::ostream &out, int DoFT, int DofCond)
2800 {
2801 
2802  int64_t nel = l2mesh->NElements();
2803  int dim = l2mesh->Dimension();
2804  TPZManVector<STATE,10> globalerrors(10,0.);
2805  for (int64_t el=0; el<nel; el++) {
2806  TPZCompEl *cel = l2mesh->ElementVec()[el];
2807  if (!cel) {
2808  continue;
2809  }
2810  TPZGeoEl *gel = cel->Reference();
2811  if (!gel || gel->Dimension() != dim) {
2812  continue;
2813  }
2814  TPZManVector<REAL,10> elerror(10,0.);
2815  elerror.Fill(0.);
2816  cel->EvaluateError(SolExataH1, elerror, 0);
2817 
2818  int nerr = elerror.size();
2819  globalerrors.resize(nerr);
2820 
2821  for (int i=0; i<nerr; i++) {
2822  globalerrors[i] += elerror[i]*elerror[i];
2823  }
2824  }
2825 
2826  out << ndiv << setw(10) << DoFT << setw(20) << DofCond << setw(28) << sqrt(globalerrors[1]) << setw(35) << sqrt(globalerrors[2]) << endl;
2827 
2828 }
2829 
2830 
2831 void hdivCurvedJCompAppMath::ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
2832 {
2833  int64_t nel = hdivmesh->NElements();
2834  int dim = hdivmesh->Dimension();
2835  TPZManVector<STATE,10> globalerrorsDual(10,0.);
2836  for (int64_t el=0; el<nel; el++) {
2837  TPZCompEl *cel = hdivmesh->ElementVec()[el];
2838  if(cel->Reference()->Dimension()!=dim) continue;
2839  TPZManVector<REAL,10> elerror(10,0.);
2840  elerror.Fill(0.);
2841  cel->EvaluateError(SolExata, elerror, 0);
2842  int nerr = elerror.size();
2843  for (int i=0; i<nerr; i++) {
2844  globalerrorsDual[i] += elerror[i]*elerror[i];
2845  }
2846 
2847 
2848  }
2849 
2850 
2851  nel = l2mesh->NElements();
2852 
2853  TPZManVector<STATE,10> globalerrorsPrimal(10,0.);
2854  for (int64_t el=0; el<nel; el++) {
2855  TPZCompEl *cel = l2mesh->ElementVec()[el];
2856  TPZManVector<REAL,10> elerror(10,0.);
2857  cel->EvaluateError(SolExata, elerror, 0);
2858  int nerr = elerror.size();
2859  globalerrorsPrimal.resize(nerr);
2860 
2861  for (int i=0; i<nerr; i++) {
2862  globalerrorsPrimal[i] += elerror[i]*elerror[i];
2863  }
2864 
2865  }
2866 
2867  // sequence of storage
2868  // for each
2869  errors.PutVal(pos, 0, p); // polinomial order
2870  errors.PutVal(pos, 1, ndiv); // refinement
2871  errors.PutVal(pos, 2, sqrt(globalerrorsPrimal[1])); // pressure
2872  errors.PutVal(pos, 3, sqrt(globalerrorsDual[1])); // flux
2873 
2874 }
2875 
2876 void hdivCurvedJCompAppMath::ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, std::ostream &out, int DoFT, int DofCond)
2877 {
2878  int64_t nel = hdivmesh->NElements();
2879  int dim = hdivmesh->Dimension();
2880  TPZManVector<STATE,10> globalerrorsDual(10,0.);
2881  for (int64_t el=0; el<nel; el++) {
2882  TPZCompEl *cel = hdivmesh->ElementVec()[el];
2883  if(cel->Reference()->Dimension()!=dim) continue;
2884  TPZManVector<REAL,10> elerror(10,0.);
2885  elerror.Fill(0.);
2886  cel->EvaluateError(SolExata, elerror, 0);
2887  int nerr = elerror.size();
2888  for (int i=0; i<nerr; i++) {
2889  globalerrorsDual[i] += elerror[i]*elerror[i];
2890  }
2891 
2892 
2893  }
2894 
2895 
2896  nel = l2mesh->NElements();
2897 
2898  TPZManVector<STATE,10> globalerrorsPrimal(10,0.);
2899  for (int64_t el=0; el<nel; el++) {
2900  TPZCompEl *cel = l2mesh->ElementVec()[el];
2901  TPZManVector<REAL,10> elerror(10,0.);
2902  cel->EvaluateError(SolExata, elerror, 0);
2903  int nerr = elerror.size();
2904  globalerrorsPrimal.resize(nerr);
2905 
2906  for (int i=0; i<nerr; i++) {
2907  globalerrorsPrimal[i] += elerror[i]*elerror[i];
2908  }
2909 
2910  }
2911 
2912  out << ndiv << setw(10) << DoFT << setw(20) << DofCond << setw(28) << sqrt(globalerrorsPrimal[1]) << setw(35) << sqrt(globalerrorsDual[1]) << endl;
2913 
2914 }
2915 
2916 void hdivCurvedJCompAppMath::ChangeExternalOrderConnects(TPZCompMesh *mesh){
2917 
2918  int nEl= mesh-> NElements();
2919  int dim = mesh->Dimension();
2920 
2921  int cordermin = -1;
2922  for (int iel=0; iel<nEl; iel++) {
2923  TPZCompEl *cel = mesh->ElementVec()[iel];
2924  if (!cel) continue;
2925  int ncon = cel->NConnects();
2926  int corder = 0;
2927  int nshape = 0;
2928 
2929  if(cel->Dimension()== dim){
2930  for (int icon=0; icon<ncon-1; icon++){
2931  TPZConnect &co = cel->Connect(icon);
2932  corder = co.Order();
2933  nshape = co.NShape();
2934  if(corder!=cordermin){
2935  cordermin = corder-1;
2936  co.SetOrder(cordermin,cel->ConnectIndex(icon));
2937  co.SetNShape(nshape-1);
2938  mesh->Block().Set(co.SequenceNumber(),nshape-1);
2939  }
2940  }
2941  }
2942  }
2943  mesh->ExpandSolution();
2944  mesh->CleanUpUnconnectedNodes();
2945 }
2946 
2947 
2948 
2949 void hdivCurvedJCompAppMath::SolveSyst(TPZAnalysis &an, TPZCompMesh *fCmesh)
2950 {
2951  std::cout <<"Numero de equacoes "<< fCmesh->NEquations()<< std::endl;
2952 
2953  bool isdirect = true;
2954  bool simetrico = true;
2955  bool isfrontal = false;
2956  if (isdirect)
2957  {
2958  if (simetrico)
2959  {
2960  //TPZSkylineStructMatrix strmat(fCmesh);
2961  if (isfrontal) {
2962  TPZParFrontStructMatrix<TPZFrontSym<STATE> > strmat(fCmesh);
2963  strmat.SetDecomposeType(ELDLt);
2964 
2965  int numthreads = 1;
2966 
2967  strmat.SetNumThreads(numthreads);
2968 
2969  an.SetStructuralMatrix(strmat);
2970  }
2971  else
2972  {
2973  //TPZBandStructMatrix full(fCmesh);
2974  TPZSkylineStructMatrix skylstr(fCmesh); //caso simetrico
2975  // TPZSkylineNSymStructMatrix full(fCmesh);
2976  an.SetStructuralMatrix(skylstr);
2977  }
2978 
2979 
2980  TPZStepSolver<STATE> step;
2981  step.SetDirect(ELDLt); //caso simetrico
2982  an.SetSolver(step);
2983  an.Run();
2984  }
2985  else
2986  {
2987  TPZBandStructMatrix full(fCmesh);
2988  an.SetStructuralMatrix(full);
2989  TPZStepSolver<STATE> step;
2990  step.SetDirect(ELU);
2991  an.SetSolver(step);
2992  an.Run();
2993  }
2994 
2995  }
2996  else
2997  {
2998  TPZSkylineStructMatrix skylstr(fCmesh); //caso simetrico
2999  skylstr.SetNumThreads(10);
3000  an.SetStructuralMatrix(skylstr);
3001 
3002  TPZAutoPointer<TPZMatrix<STATE> > matbeingcopied = skylstr.Create();
3003  TPZAutoPointer<TPZMatrix<STATE> > matClone = matbeingcopied->Clone();
3004 
3005  TPZStepSolver<STATE> *precond = new TPZStepSolver<STATE>(matClone);
3006  TPZStepSolver<STATE> *Solver = new TPZStepSolver<STATE>(matbeingcopied);
3007  precond->SetReferenceMatrix(matbeingcopied);
3008  precond->SetDirect(ELDLt);
3009  Solver->SetGMRES(20, 20, *precond, 1.e-18, 0);
3010  // Solver->SetCG(10, *precond, 1.0e-10, 0);
3011  an.SetSolver(*Solver);
3012  an.Run();
3013  }
3014 
3015 
3016 
3017 }
3018 
3019 
3020 
3021 
3022 
3023 
RotationMatrix
void RotationMatrix(TPZFMatrix< REAL > &R, REAL thetaRad, int axis)
Definition: TPZPlasticityTest.h:2349
TPZCompMesh::NElements
int64_t NElements() const
Number of computational elements allocated.
Definition: pzcmesh.h:169
TPZCompMesh::Reference
TPZGeoMesh * Reference() const
Returns a pointer to the geometrical mesh associated.
Definition: pzcmesh.h:209
hdivCurvedJCompAppMath::ForcingBC2D
static void ForcingBC2D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdivCurvedJCompAppMath.cpp:1511
atan2
TPZFlopCounter atan2(const TPZFlopCounter &val1, const TPZFlopCounter &val2)
Returns the arc tangent in radians and increments the counter of the Arc Tangent. ATAN2 returns the ...
Definition: pzreal.h:544
TPZLagrangeMultiplier
Material which implements a Lagrange Multiplier.
Definition: TPZLagrangeMultiplier.h:16
TPZMatMixedPoisson3D::SetPermeabilityTensor
void SetPermeabilityTensor(TPZFMatrix< REAL > K, TPZFMatrix< REAL > invK)
Definition: pzmatmixedpoisson3d.h:115
TPZCompMesh::AdjustBoundaryElements
void AdjustBoundaryElements()
Will refine the elements associated with a boundary condition till there are no elements constrained ...
Definition: pzcmesh.cpp:1423
TPZGeoNode::SetCoord
void SetCoord(const TPZVec< REAL > &x)
Sets all coordinates into the current node. It gets the dim values from x.
Definition: pzgnode.cpp:60
TPZCompMesh::SetAllCreateFunctionsContinuous
void SetAllCreateFunctionsContinuous()
Definition: pzcmesh.h:508
hdivCurvedJCompAppMath::ParametricSphere
TPZManVector< REAL, 3 > ParametricSphere(REAL radius, REAL phi, REAL theta)
Definition: hdivCurvedJCompAppMath.cpp:926
TPZConnect::IncrementElConnected
void IncrementElConnected()
Increment fNElConnected.
Definition: pzconnect.h:260
TPZCompElSide
Implements computational element and a side. Computational Element.
Definition: pzcompel.h:632
TPZBlock::Set
int Set(const int index, const int dim, const int pos=-1)
Modifies existing block dimensions or creates a new block with given index.
Definition: pzblock.cpp:104
TPZConnect
Represents a set of shape functions associated with a computational element/side. Computational Eleme...
Definition: pzconnect.h:30
ELDLt
Definition: pzmatrix.h:52
TPZAnalysis::SetStructuralMatrix
void SetStructuralMatrix(TPZAutoPointer< TPZStructMatrix > strmatrix)
Set structural matrix as auto pointer for analysis.
Definition: pzanalysis.cpp:82
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
TPZBandStructMatrix
Implements Banded Structural Matrices. Structural Matrix.
Definition: pzbstrmatrix.h:18
hdivCurvedJCompAppMath::CMeshFlux
TPZCompMesh * CMeshFlux(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdivCurvedJCompAppMath.cpp:1757
TPZCompEl::SetConnectIndex
virtual void SetConnectIndex(int inode, int64_t index)=0
Set the index i to node inode.
TPZGeoEl::MaterialId
int MaterialId() const
Returns the material index of the element.
Definition: pzgeoel.h:250
TPZCompEl::EvaluateError
virtual void EvaluateError(std::function< void(const TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv)> func, TPZVec< REAL > &errors, bool store_error)
Performs an error estimate on the elemen.
Definition: pzcompel.cpp:415
TPZElementGroup::NConnects
virtual int NConnects() const override
Returns the number of nodes of the element.
Definition: pzelementgroup.h:93
TPZVec::resize
virtual void resize(const int64_t newsize)
Definition: pzvec.h:213
TPZBuildMultiphysicsMesh::AddElements
static void AddElements(TPZVec< TPZCompMesh *> &cmeshVec, TPZCompMesh *MFMesh)
Creating multiphysic elements into mphysics computational mesh. Method to add elements in the mesh mu...
Definition: pzbuildmultiphysicsmesh.cpp:40
TPZConnect::SetOrder
void SetOrder(int order, int64_t index)
Set the order of the shapefunction associated with the connect.
Definition: pzconnect.h:168
TPZCompMesh::NEquations
int64_t NEquations()
This computes the number of equations associated with non-restrained nodes.
Definition: pzcmesh.cpp:721
TPZCompMesh::AllocateNewConnect
virtual int64_t AllocateNewConnect(int nshape, int nstate, int order)
Returns an index to a new connect.
Definition: pzcmesh.h:693
TPZStepSolver
Defines step solvers class. Solver.
Definition: pzmganalysis.h:17
TPZAnalysis::SetSolver
void SetSolver(TPZMatrixSolver< STATE > &solver)
Set solver matrix.
Definition: pzanalysis.cpp:1202
TPZCompEl::Dimension
virtual int Dimension() const =0
Dimension of the element.
hdivCurvedJCompAppMath::MakeCircle
TPZGeoMesh * MakeCircle(int ndiv)
Definition: hdivCurvedJCompAppMath.cpp:491
TPZBuildMultiphysicsMesh::TransferFromMultiPhysics
static void TransferFromMultiPhysics(TPZVec< TPZCompMesh *> &cmeshVec, TPZCompMesh *MFMesh)
Transfer information from a specific mesh multiphysics for the current specific set of meshes...
Definition: pzbuildmultiphysicsmesh.cpp:400
TPZConnect::SetLagrangeMultiplier
void SetLagrangeMultiplier(unsigned char mult)
Set the connect as a pressure connect or not.
Definition: pzconnect.h:236
TPZGeoMesh::NElements
int64_t NElements() const
Number of elements of the mesh.
Definition: pzgmesh.h:129
TPZGeoElSide
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
TPZElementGroup
Class which groups elements to characterize dense matrices.
Definition: pzelementgroup.h:20
hdivCurvedJCompAppMath::ErrorPrimalDual
void ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
Definition: hdivCurvedJCompAppMath.cpp:2831
hdivCurvedJCompAppMath::SolExataH1
static void SolExataH1(const TPZVec< REAL > &pt, TPZVec< STATE > &solp, TPZFMatrix< STATE > &flux)
Definition: hdivCurvedJCompAppMath.cpp:1402
hdivCurvedJCompAppMath::ErrorH1
void ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
Definition: hdivCurvedJCompAppMath.cpp:2764
TPZChunkVector::NElements
int64_t NElements() const
Access method to query the number of elements of the vector.
Definition: TPZChunkVector.h:236
TPZGeoElSide::Neighbour
TPZGeoElSide Neighbour() const
Definition: pzgeoel.h:754
TPZGeoEl::NSides
virtual int NSides() const =0
Returns the number of connectivities of the element.
TPZGeoEl::Divide
virtual void Divide(TPZVec< TPZGeoEl *> &pv)
Divides the element and puts the resulting elements in the vector.
Definition: pzgeoel.h:742
hdivCurvedJCompAppMath::ForcingBC4D
static void ForcingBC4D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdivCurvedJCompAppMath.cpp:1617
TPZGeoMesh::SetDimension
void SetDimension(int dim)
Set Dimension.
Definition: pzgmesh.h:285
norder
const int norder
Definition: hdivCurvedJCompAppMath.cpp:11
TPZCompMesh::LoadReferences
void LoadReferences()
Map this grid in the geometric grid.
Definition: pzcmesh.cpp:482
TPZCreateApproximationSpace::CreateDisconnectedElements
void CreateDisconnectedElements(bool create)
Determine if the mesh will be created with disconnected elements After the mesh is created...
Definition: pzcreateapproxspace.h:165
TPZBuildMultiphysicsMesh::AddConnects
static void AddConnects(TPZVec< TPZCompMesh *> &cmeshVec, TPZCompMesh *MFMesh)
Definition: pzbuildmultiphysicsmesh.cpp:108
TPZCompMesh::SetDefaultOrder
void SetDefaultOrder(int order)
Definition: pzcmesh.h:403
TPZSkylineStructMatrix
Implements SkyLine Structural Matrices. Structural Matrix.
Definition: pzskylstrmatrix.h:19
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
hdivCurvedJCompAppMath::PrintErrors
void PrintErrors(geomDomain geodomain, ApproximationSpace problem, Eltype element, TPZVec< int > POrderBeginAndEnd, TPZVec< int > ndivinterval, TPZVec< REAL > &errors, std::ostream &output)
Definition: hdivCurvedJCompAppMath.cpp:299
TPZConnect::LagrangeMultiplier
unsigned char LagrangeMultiplier() const
Access method to return the indication whether the connect is associated with a pressure lagrange mul...
Definition: pzconnect.h:230
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZMaterial::SetForcingFunctionExact
void SetForcingFunctionExact(TPZAutoPointer< TPZFunction< STATE > > fp)
Sets a procedure as exact solution for the problem.
Definition: TPZMaterial.h:405
acos
expr_ expr_ expr_ expr_ acos
Definition: tfadfunc.h:75
TPZCompMesh::Element
TPZCompEl * Element(int64_t iel)
Definition: pzcmesh.h:185
TPZAdmChunkVector::Resize
void Resize(const int newsize)
Increase the size of the chunk vector.
Definition: pzadmchunk.h:280
sin
sin
Definition: fadfunc.h:63
TPZConnect::Order
unsigned char Order() const
Access function to return the order associated with the connect.
Definition: pzconnect.h:217
TPZCompMesh::ApproxSpace
TPZCreateApproximationSpace & ApproxSpace()
Definition: pzcmesh.h:498
TPZConnect::SequenceNumber
int64_t SequenceNumber() const
Returns the Sequence number of the connect object.
Definition: pzconnect.h:158
TPZCompElDisc::SetTensorialShape
static void SetTensorialShape(TPZCompMesh *cmesh)
Sets tensorial shape functions for all Discontinuous elements in cmesh.
Definition: TPZCompElDisc.cpp:77
TPZAnalysis
Implements the sequence of actions to perform a finite element analysis. Analysis.
Definition: pzanalysis.h:32
TPZSkylineStructMatrix::Create
virtual TPZMatrix< STATE > * Create()
Definition: pzskylstrmatrix.cpp:33
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
TPZMatMixedPoisson3D::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: pzmatmixedpoisson3d.cpp:114
TPZFrontStructMatrix::SetDecomposeType
virtual void SetDecomposeType(DecomposeType dectype)
Set the decomposition type.
Definition: TPZFrontStructMatrix.h:67
TPZGeoMesh::SetMaxNodeId
void SetMaxNodeId(int64_t id)
Used in patch meshes.
Definition: pzgmesh.h:120
hdivCurvedJCompAppMath::ForcingH1
static void ForcingH1(const TPZVec< REAL > &pt, TPZVec< STATE > &ff, TPZFMatrix< STATE > &flux)
Definition: hdivCurvedJCompAppMath.cpp:1432
TPZGeoElRefPattern
Implements a generic geometric element which is refined according to a generic refinement pattern...
Definition: pzgmesh.h:35
TPZBuildMultiphysicsMesh::TransferFromMeshes
static void TransferFromMeshes(TPZVec< TPZCompMesh *> &cmeshVec, TPZCompMesh *MFMesh)
Transfer information from a specific set of meshes for the current mesh multiphysics.
Definition: pzbuildmultiphysicsmesh.cpp:321
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
TPZInterpolationSpace::NSideConnects
virtual int NSideConnects(int iside) const =0
Returns the number of dof nodes along side iside.
TPZCompMesh::SetAllCreateFunctionsDiscontinuous
void SetAllCreateFunctionsDiscontinuous()
Definition: pzcmesh.h:503
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
hdivCurvedJCompAppMath::RotateGeomesh
void RotateGeomesh(TPZGeoMesh *gmesh, REAL CounterClockwiseAngle, int &Axis)
Definition: hdivCurvedJCompAppMath.cpp:859
TPZGeoMesh::NNodes
int64_t NNodes() const
Number of nodes of the mesh.
Definition: pzgmesh.h:126
TPZCompMesh::ExpandSolution
virtual void ExpandSolution()
Adapt the solution vector to new block dimensions.
Definition: pzcmesh.cpp:396
TPZStructMatrixBase::SetNumThreads
virtual void SetNumThreads(int n)
Definition: TPZStructMatrixBase.h:41
hdivCurvedJCompAppMath::CMeshMixed
TPZCompMesh * CMeshMixed(TPZGeoMesh *gmesh, TPZVec< TPZCompMesh *> meshvec)
Definition: hdivCurvedJCompAppMath.cpp:2189
TPZCompMesh::SetAllCreateFunctionsHDiv
void SetAllCreateFunctionsHDiv()
Definition: pzcmesh.h:523
substruct_tst14.test.numthreads
list numthreads
Definition: test.py:210
ELU
Definition: pzmatrix.h:52
TPZCompMesh::Dimension
int Dimension() const
Returns the dimension of the simulation.
Definition: pzcmesh.h:148
TPZCompMesh::ComputeNodElCon
virtual void ComputeNodElCon()
Compute the number of elements connected to each connect object.
Definition: pzcmesh.cpp:668
TPZGeoMesh::NodeVec
TPZAdmChunkVector< TPZGeoNode > & NodeVec()
Definition: pzgmesh.h:140
TPZGeoEl::Reference
TPZCompEl * Reference() const
Return a pointer to the element referenced by the geometric element.
Definition: pzgeoel.h:750
TPZCompMesh::AutoBuild
virtual void AutoBuild(const std::set< int > *MaterialIDs)
Creates the computational elements, and the degree of freedom nodes.
Definition: pzcmesh.cpp:308
TPZConnect::NShape
unsigned int NShape() const
Definition: pzconnect.h:151
hdivCurvedJCompAppMath::ParametricCircle
TPZManVector< REAL, 3 > ParametricCircle(REAL radius, REAL theta)
Definition: hdivCurvedJCompAppMath.cpp:682
TPZConnect::SetNShape
void SetNShape(int nshape)
Set the number of shape functions associated with the connect.
Definition: pzconnect.h:205
TPZMatrixSolver::SetReferenceMatrix
virtual void SetReferenceMatrix(TPZAutoPointer< TPZMatrix< TVar > > matrix)
This method gives a preconditioner to share a matrix with the referring solver object.
Definition: pzsolve.h:132
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
hdivCurvedJCompAppMath::CMeshH1
TPZCompMesh * CMeshH1(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdivCurvedJCompAppMath.cpp:1700
hdivCurvedJCompAppMath::SolveSyst
void SolveSyst(TPZAnalysis &an, TPZCompMesh *fCmesh)
Definition: hdivCurvedJCompAppMath.cpp:2949
TPZCompMesh::Block
const TPZBlock< STATE > & Block() const
Access the block structure of the solution vector.
Definition: pzcmesh.h:213
TPZCompMesh::NConnects
int64_t NConnects() const
Number of connects allocated including free nodes.
Definition: pzcmesh.h:163
co
REAL co[8][3]
Coordinates of the eight nodes.
Definition: tpzgensubstruct.cpp:56
TPZCompMesh::ConnectVec
TPZAdmChunkVector< TPZConnect > & ConnectVec()
Return a reference to the connect pointers vector.
Definition: pzcmesh.h:198
TPZCondensedCompEl
Class which implements an element which condenses the internal connects.
Definition: pzcondensedcompel.h:32
TPZCompEl::Index
int64_t Index() const
Returns element index of the mesh fELementVec list.
Definition: pzcompel.h:821
TPZVTKGeoMesh::PrintGMeshVTK
static void PrintGMeshVTK(TPZGeoMesh *gmesh, std::ofstream &file, bool matColor=true)
Default constructor for graphical mesh with VTK format.
Definition: TPZVTKGeoMesh.cpp:197
hdivCurvedJCompAppMath::Forcing
static void Forcing(const TPZVec< REAL > &pt, TPZVec< STATE > &ff)
Definition: hdivCurvedJCompAppMath.cpp:1350
TPZCompEl::Connect
virtual TPZConnect & Connect(int i) const
Returns a pointer to the ith node.
Definition: pzcompel.cpp:298
TPZCompMesh::SetDimModel
void SetDimModel(int dim)
Set de dimension of the domain of the problem.
Definition: pzcmesh.h:139
TPZAnalysis::Run
virtual void Run(std::ostream &out=std::cout)
Calls the appropriate sequence of methods to build a solution or a time stepping sequence.
Definition: pzanalysis.cpp:920
TPZCompEl::ConnectIndex
virtual int64_t ConnectIndex(int i) const =0
Returns the index of the ith connectivity of the element.
TPZInterfaceElement
Computes the contribution over an interface between two discontinuous elements. Computational Element...
Definition: TPZInterfaceEl.h:30
TPZGeoEl::HasSubElement
virtual int HasSubElement() const =0
Return 1 if the element has subelements.
hdivCurvedJCompAppMath::ForcingBC3D
static void ForcingBC3D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdivCurvedJCompAppMath.cpp:1564
hdivCurvedJCompAppMath::Run
void Run(geomDomain geodomain, ApproximationSpace problem, Eltype element, TPZVec< int > POrderBeginAndEnd, TPZVec< int > ndivinterval, TPZFMatrix< REAL > &errors)
Definition: hdivCurvedJCompAppMath.cpp:107
angle
static REAL angle
Angle in radians to test.
Definition: pzsubcmesh.cpp:53
TPZInterpolationSpace
Implements the interfaces for TPZCompElDisc, TPZInterfaceElement and TPZInterpolatedElement. Computational element.
Definition: pzinterpolationspace.h:17
hdivCurvedJCompAppMath::SetupDisconnectedHdivboud
void SetupDisconnectedHdivboud(const int left, const int rigth, TPZCompMesh *cmesh)
Definition: hdivCurvedJCompAppMath.cpp:1895
hdivCurvedJCompAppMath::GMeshCilindricalMesh
TPZGeoMesh * GMeshCilindricalMesh(int ndiv)
Definition: hdivCurvedJCompAppMath.cpp:935
TPZCompEl::NConnects
virtual int NConnects() const =0
Returns the number of nodes of the element.
hdivCurvedJCompAppMath::ForcingBC1D
static void ForcingBC1D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdivCurvedJCompAppMath.cpp:1460
TPZGeoElSide::Element
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
hdivCurvedJCompAppMath::MakeSphereFromQuadrilateral
TPZGeoMesh * MakeSphereFromQuadrilateral(int dimensao, bool triang, int ndiv)
Definition: hdivCurvedJCompAppMath.cpp:691
TPZGeoMesh::BuildConnectivity
void BuildConnectivity()
Build the connectivity of the grid.
Definition: pzgmesh.cpp:967
TPZMaterial::CreateBC
virtual TPZBndCond * CreateBC(TPZMaterial *reference, int id, int typ, TPZFMatrix< STATE > &val1, TPZFMatrix< STATE > &val2)
Creates an object TPZBndCond derived of TPZMaterial.
Definition: TPZMaterial.cpp:282
TPZInterpolationSpace::SetSideOrient
virtual void SetSideOrient(int side, int sideorient)
It set the normal orientation of the element by the side. Only side that has dimension equal to my di...
Definition: pzinterpolationspace.cpp:1718
TPZCompMesh::InsertMaterialObject
int InsertMaterialObject(TPZMaterial *mat)
Insert a material object in the datastructure.
Definition: pzcmesh.cpp:287
TPZCompElDisc::SetTotalOrderShape
static void SetTotalOrderShape(TPZCompMesh *cmesh)
Set total order shape functions for all Discontinuous elements in cmesh.
Definition: TPZCompElDisc.cpp:89
TPZInterpolationSpace::SideConnectLocId
virtual int SideConnectLocId(int icon, int is) const =0
Returns the local node number of icon along is.
TPZMatPoisson3d
DESCRIBE PLEASE.
Definition: pzpoisson3d.h:26
TPZCompMesh::CleanUpUnconnectedNodes
virtual void CleanUpUnconnectedNodes()
Delete the nodes which have no elements connected to them.
Definition: pzcmesh.cpp:498
TPZGeoEl::Dimension
virtual int Dimension() const =0
Returns the dimension of the element.
TPZCompElSide::Element
TPZCompEl * Element() const
Gives a pointer to the reference computational element.
Definition: pzcompel.h:675
TPZGeoMesh::Dimension
int Dimension()
Get Dimension.
Definition: pzgmesh.h:291
TPZCompEl::Reference
TPZGeoEl * Reference() const
Return a pointer to the corresponding geometric element if such exists, return 0 otherwise.
Definition: pzcompel.cpp:1137
TPZGeoNode
Implements a geometric node in the pz environment. Geometry.
Definition: pzgnode.h:31
atan
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ atan
Definition: tfadfunc.h:85
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
TPZMaterial::SetForcingFunction
void SetForcingFunction(TPZAutoPointer< TPZFunction< STATE > > fp)
Sets a procedure as source function for the material.
Definition: TPZMaterial.h:368
TPZDummyFunction::SetPolynomialOrder
void SetPolynomialOrder(int porder)
Definition: pzfunction.h:226
TPZCompMesh
Implements computational mesh. Computational Mesh.
Definition: pzcmesh.h:47
TPZCompMesh::ElementVec
TPZAdmChunkVector< TPZCompEl * > & ElementVec()
Returns a reference to the element pointers vector.
Definition: pzcmesh.h:183
TPZCompElDisc::SetConstC
void SetConstC(REAL c)
Definition: TPZCompElDisc.h:216
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
TPZCompElSide::Side
int Side() const
Returns the side index.
Definition: pzcompel.h:688
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
hdivCurvedJCompAppMath::SolExata
static void SolExata(const TPZVec< REAL > &pt, TPZVec< STATE > &solp, TPZFMatrix< STATE > &flux)
Definition: hdivCurvedJCompAppMath.cpp:1226
hdivCurvedJCompAppMath::ChangeExternalOrderConnects
void ChangeExternalOrderConnects(TPZCompMesh *mesh)
Definition: hdivCurvedJCompAppMath.cpp:2916
TPZCompElDisc::SetCenterPoint
void SetCenterPoint(int i, REAL x)
Definition: TPZCompElDisc.h:314
TPZBuildMultiphysicsMesh::AddWrap
static void AddWrap(TPZMultiphysicsElement *mfcel, int matskeleton, TPZStack< TPZStack< TPZMultiphysicsElement *, 7 > > &ListGroupEl)
Create skeleton elements of the wrap of me.
Definition: pzbuildmultiphysicsmesh.cpp:798
TPZMatPoisson3d::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: pzpoisson3d.h:161
TPZStepSolver::SetDirect
void SetDirect(const DecomposeType decomp)
Definition: pzstepsolver.cpp:184
TPZCompElDisc::SetTrueUseQsiEta
void SetTrueUseQsiEta()
Definition: TPZCompElDisc.h:218
TPZCompElDisc
This class implements a discontinuous element (for use with discontinuous Galerkin). Computational Element.
Definition: TPZCompElDisc.h:35
TPZGeoNode::GetCoordinates
void GetCoordinates(TPZVec< REAL > &co)
Fill the coordinates of the node.
Definition: pzgnode.cpp:83
hdivCurvedJCompAppMath::ForcingBC5D
static void ForcingBC5D(const TPZVec< REAL > &pt, TPZVec< STATE > &solp)
Definition: hdivCurvedJCompAppMath.cpp:1666
TPZGeoMesh::ResetReference
void ResetReference()
Resets all load references in elements and nodes.
Definition: pzgmesh.cpp:197
TPZParFrontStructMatrix
Is a structural matrix with parallel techniques included. Structural Matrix Frontal.
Definition: TPZParFrontStructMatrix.h:42
hdivCurvedJCompAppMath::CMeshPressure
TPZCompMesh * CMeshPressure(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdivCurvedJCompAppMath.cpp:1957
TPZCompEl
Defines the interface of a computational element. Computational Element.
Definition: pzcompel.h:59
TPZMatMixedPoisson3D
Material to solve a mixed poisson problem 3D by multiphysics simulation.
Definition: pzmatmixedpoisson3d.h:33
TPZCompMesh::SetAllCreateFunctionsMultiphysicElem
void SetAllCreateFunctionsMultiphysicElem()
Definition: pzcmesh.h:542
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZFMatrix::PutVal
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
TPZElementGroup::AddElement
virtual void AddElement(TPZCompEl *cel)
add an element to the element group
Definition: pzelementgroup.cpp:45
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZStepSolver::SetGMRES
void SetGMRES(const int64_t numiterations, const int numvectors, const TPZMatrixSolver< TVar > &pre, const REAL tol, const int64_t FromCurrent)
Definition: pzstepsolver.cpp:202
TPZGeoMesh::ElementVec
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138
hdivCurvedJCompAppMath::RotateNode
void RotateNode(TPZVec< REAL > &iCoords, REAL CounterClockwiseAngle, int &Axis)
Definition: hdivCurvedJCompAppMath.cpp:1172
TPZAutoPointer
This class implements a reference counter mechanism to administer a dynamically allocated object...
Definition: TPZPersistenceManager.h:14
Generated on Sun Aug 16 2020 11:17:09 for NeoPZ by   doxygen 1.8.13