NeoPZ
hdiv3dpaper201504.cpp
Go to the documentation of this file.
1 //
2 // hdiv3dpaper201504.cpp
3 // PZ
4 //
5 // Created by Douglas Castro on 22/04/15.
6 //
7 //
8 
9 #include "hdiv3dpaper201504.h"
10 
11 
12 Hdiv3dPaper201504::Hdiv3dPaper201504()
13 {
14 
15  fDim = 3;
16 
17  fmatId = 1;
18 
19  fdirichlet = 0;
20  fneumann = 1;
21 
22  fbc0 = -1;
23  fbc1 = -2;
24  fbc2 = -3;
25  fbc3 = -4;
26  fbc4 = -5;
27  fbc5 = -6;
28  fmatskeleton = -7;
29 
30  fisH1 = false;
31 
32  ftetra = false;
33 
34  fprisma = false;
35 
36  tetraedra_2.Resize(6, 4);
37 
38  tetraedra_2(0,0) = 1;
39  tetraedra_2(0,1) = 2;
40  tetraedra_2(0,2) = 5;
41  tetraedra_2(0,3) = 4;
42 
43  tetraedra_2(1,0) = 4;
44  tetraedra_2(1,1) = 7;
45  tetraedra_2(1,2) = 3;
46  tetraedra_2(1,3) = 2;
47 
48  tetraedra_2(2,0) = 0;
49  tetraedra_2(2,1) = 1;
50  tetraedra_2(2,2) = 2;
51  tetraedra_2(2,3) = 4;
52 
53  tetraedra_2(3,0) = 0;
54  tetraedra_2(3,1) = 2;
55  tetraedra_2(3,2) = 3;
56  tetraedra_2(3,3) = 4;
57 
58  tetraedra_2(4,0) = 4;
59  tetraedra_2(4,1) = 5;
60  tetraedra_2(4,2) = 6;
61  tetraedra_2(4,3) = 2;
62 
63  tetraedra_2(5,0) = 4;
64  tetraedra_2(5,1) = 6;
65  tetraedra_2(5,2) = 7;
66  tetraedra_2(5,3) = 2;
67 
68 }
69 
70 Hdiv3dPaper201504::~Hdiv3dPaper201504()
71 {
72 
73 }
74 
75 void Hdiv3dPaper201504::Run(ApproximationSpace problem, Eltype element, TPZVec<int> POrderBeginAndEnd, TPZVec<int> ndivinterval, TPZFMatrix< REAL > &errors)
76 {
77  if (POrderBeginAndEnd.size()!=2)
78  {
79  std::cout << " POrderBeginAndEnd Vector must have size == 2 and contain just two values, the first and last approximation orders. " << std::endl;
80  std::cout << " Some like POrderBeginAndEnd[0] = 1, POrderBeginAndEnd[1] = 3 " << std::endl;
81  return;
82  }
83 
84  if (ndivinterval.size()!=2)
85  {
86  std::cout << " ndivinterval Vector must have size == 2 and contain just two values, the first and last number of uniform refinements. " << std::endl;
87  std::cout << " Some like ndivinterval[0] = 0, ndivinterval[1] = 3 for simulations with 0, 1, 2 and 3 refinements " << std::endl;
88  std::cout << " Some like ndivinterval[0] = 1, ndivinterval[1] = 1 for simulations with 1 refinement " << std::endl;
89  return;
90  }
91 
92  int sizeErrorMatrix = (ndivinterval[1]-ndivinterval[0])*(POrderBeginAndEnd[1] - POrderBeginAndEnd[0]) + 1;
93  int errorposition = 0;
94 
95  errors.Resize(sizeErrorMatrix, 4);
96 
97  for(int p = POrderBeginAndEnd[0]; p <= POrderBeginAndEnd[1]; p++)
98  {
99 
100  for (int ndiv = ndivinterval[0]; ndiv<=ndivinterval[1]; ndiv++)
101  {
102  int ordemP = p;
103 
104  std::cout<< " BEGIN (CUBE Domain) - Polinomial degree: " << ordemP << " rerfinement size (h): " << ndiv << std::endl;
105 
106  TPZGeoMesh *gmesh;
107  switch (element) {
108  case ECub:
109  {
110  gmesh = this->CreateOneCubo(ndiv);
111  //gmesh = this->CreateOneQuadraticCube(ndiv);
112  }
113  break;
114  case EPrism:
115  {
116  gmesh = this->GMeshWithPrism( ndiv);
117  }
118  break;
119  case ETetra:
120  {
121  double dndiv = ndiv;
122  int nref = (int) pow(2., dndiv);
123  gmesh = this->CreateOneCuboWithTetraedrons(nref);
124  }
125  break;
126  default:
127  break;
128  }
129 
130 
131 
132  switch (problem) {
133  case EH1:
134  {
135 
136  gmesh->SetDimension(fDim);
137  TPZCompMesh *cmeshH1 = this->CMeshH1(gmesh, ordemP, fDim);
138  //condensar
139  for (int64_t iel=0; iel<cmeshH1->NElements(); iel++) {
140  TPZCompEl *cel = cmeshH1->Element(iel);
141  if(!cel) continue;
142  new TPZCondensedCompEl(cel);
143  }
144 
145  cmeshH1->ExpandSolution();
146  cmeshH1->CleanUpUnconnectedNodes();
147 
148  TPZAnalysis anh1(cmeshH1, true);
149 
150  SolveSyst(anh1, cmeshH1);
151 
152  ErrorH1(cmeshH1, ordemP, ndiv, errorposition, errors);
153 
154  }
155  break;
156  case EHDiv:
157  {
158  DebugStop(); //Mudanca na topologia
159  }
160  break;
161  case EHDivStar:
162  {
163 
164  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
165 
166  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
167 
168  // P**
169  ChangeExternalOrderConnects(cmesh1);
170 
171 
172  //malha multifisica
173  TPZVec<TPZCompMesh *> meshvec(2);
174  meshvec[0] = cmesh1;
175  meshvec[1] = cmesh2;
176 
177  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
178 
179  TPZAnalysis an(mphysics, true);
180 
181  SolveSyst(an, mphysics);
182 
183  //Calculo do erro
184  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
185  TPZVec<REAL> erros;
186 
187  std::cout << "Computing Errors\n";
188 
189  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, errorposition, errors);
190  }
191  break;
192 
193  case EHDivStarStar:
194  {
195 
196  ordemP = ordemP + 1;
197 
198  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
199 
200  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
201  //P**
202  ChangeExternalOrderConnects(cmesh1);
203 
204  //malha multifisica
205  TPZVec<TPZCompMesh *> meshvec(2);
206  meshvec[0] = cmesh1;
207  meshvec[1] = cmesh2;
208 
209  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
210 
211  TPZAnalysis an(mphysics, true);
212 
213  SolveSyst(an, mphysics);
214 
215  //Calculo do erro
216  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
217  TPZVec<REAL> erros;
218 
219  std::cout << "Computing Errors\n";
220 
221  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, errorposition, errors);
222  }
223  break;
224 
225  default:
226  break;
227  }
228 
229  std::cout<< " END - polinomial degree: " << ordemP << " refinement size (h): " << ndiv << std::endl;
230 
231  errorposition++;
232 
233  }
234  }
235 
236  std::cout<< " The End " << std::endl;
237 
238 
239 }
240 
241 void Hdiv3dPaper201504::PrintErrors(ApproximationSpace problem, Eltype element, TPZVec<int> POrderBeginAndEnd, TPZVec<int> ndivinterval, TPZVec<REAL> &errors, std::ostream &output)
242 {
243  if (POrderBeginAndEnd.size()!=2)
244  {
245  std::cout << " POrderBeginAndEnd Vector must have size == 2 and contain just two values, the first and last approximation orders. " << std::endl;
246  std::cout << " Some like POrderBeginAndEnd[0] = 1, POrderBeginAndEnd[1] = 3 " << std::endl;
247  return;
248  }
249 
250  if (ndivinterval.size()!=2)
251  {
252  std::cout << " ndivinterval Vector must have size == 2 and contain just two values, the first and last number of uniform refinements. " << std::endl;
253  std::cout << " Some like ndivinterval[0] = 0, ndivinterval[1] = 3 for simulations with 0, 1, 2 and 3 refinements " << std::endl;
254  std::cout << " Some like ndivinterval[0] = 1, ndivinterval[1] = 1 for simulations with 1 refinement " << std::endl;
255  return;
256  }
257 
258  for(int p = POrderBeginAndEnd[0]; p <= POrderBeginAndEnd[1]; p++)
259  {
260  output << "\n WHEN p = " << p << " \n " << endl;
261  output << "ndiv " << setw(6) << "DoFTot" << setw(20) << "DofCond" << setw(28) << "PrimalL2Error" << setw(35) << "L2DualError" << endl;
262 
263  for (int ndiv = ndivinterval[0]; ndiv<=ndivinterval[1]; ndiv++)
264  {
265  int ordemP = p;
266 
267  std::cout<< " BEGIN (Cube Domain) - Polinomial degree: " << ordemP << " rerfinement size (h): " << ndiv << std::endl;
268 
269  TPZGeoMesh *gmesh;
270  switch (element) {
271  case ECub:
272  {
273  gmesh = this->CreateOneCubo(ndiv);
274  //gmesh = this->CreateOneQuadraticCube(ndiv);
275  }
276  break;
277  case EPrism:
278  {
279  gmesh = this->GMeshWithPrism( ndiv);
280  }
281  break;
282  case ETetra:
283  {
284  double dndiv = ndiv;
285  int nref = (int) pow(2., dndiv);
286  gmesh = this->CreateOneCuboWithTetraedrons(nref);
287  }
288  break;
289  default:
290  break;
291  }
292 
293  gmesh->SetDimension(fDim);
294 
295  switch (problem) {
296  case EH1:
297  {
298 
299  TPZCompMesh *cmeshH1 = this->CMeshH1(gmesh, ordemP, fDim);
300 
301  int dofTotal = cmeshH1->NEquations();
302 
303  //condensar
304  for (int64_t iel=0; iel<cmeshH1->NElements(); iel++) {
305  TPZCompEl *cel = cmeshH1->Element(iel);
306  if(!cel) continue;
307  new TPZCondensedCompEl(cel);
308  }
309 
310  cmeshH1->ExpandSolution();
311  cmeshH1->CleanUpUnconnectedNodes();
312 
313  int dofCondensed = cmeshH1->NEquations();
314 
315  TPZAnalysis anh1(cmeshH1, true);
316 
317  SolveSyst(anh1, cmeshH1);
318 
319  ErrorH1(cmeshH1, ordemP, ndiv, output, dofTotal, dofCondensed);
320 
321  }
322  break;
323  case EHDiv:
324  {
325  DebugStop(); //Mudanca na topologia
326  }
327  break;
328  case EHDivStar:
329  {
330 
331  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
332 
333  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
334 
335  // para rodar P**
336  ChangeExternalOrderConnects(cmesh1);
337 
338  int DofCond, DoFT;
339  DoFT = cmesh1->NEquations() + cmesh2->NEquations();
340 
341  //malha multifisica
342  TPZVec<TPZCompMesh *> meshvec(2);
343  meshvec[0] = cmesh1;
344  meshvec[1] = cmesh2;
345 
346  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
347 
348  DofCond = mphysics->NEquations();
349 
350  TPZAnalysis an(mphysics, true);
351 
352  SolveSyst(an, mphysics);
353 
354  //Calculo do erro
355  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
356  TPZVec<REAL> erros;
357 
358  std::cout << "Computing Errors\n";
359 
360  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, output, DoFT, DofCond);
361  }
362  break;
363 
364  case EHDivStarStar:
365  {
366 
367 
368  ordemP = ordemP + 1;
369 
370  TPZCompMesh *cmesh1 = this->CMeshFlux(gmesh, ordemP, fDim);
371 
372  TPZCompMesh *cmesh2 = this->CMeshPressure(gmesh, ordemP, fDim);
373 
374  // para rodar P**
375  ChangeExternalOrderConnects(cmesh1);
376 
377  int DofCond, DoFT;
378  DoFT = cmesh1->NEquations() + cmesh2->NEquations();
379 
380  //malha multifisica
381  TPZVec<TPZCompMesh *> meshvec(2);
382  meshvec[0] = cmesh1;
383  meshvec[1] = cmesh2;
384 
385  TPZCompMesh * mphysics = CMeshMixed(gmesh,meshvec);
386 
387  DofCond = mphysics->NEquations();
388 
389  TPZAnalysis an(mphysics, true);
390 
391  SolveSyst(an, mphysics);
392 
393  //Calculo do erro
394  TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(meshvec, mphysics);
395  TPZVec<REAL> erros;
396 
397  std::cout << "Computing Errors\n";
398 
399  ErrorPrimalDual( cmesh2, cmesh1, ordemP, ndiv, output, DoFT, DofCond);
400  }
401  break;
402 
403  default:
404  break;
405  }
406 
407 
408  }
409 
410  output << "\n ----------------------------------------------------------------------------- " << endl;
411  std::cout<< " END (Quadrilateral Domain) - Polinomial degree: " << p << std::endl;
412 
413  }
414 
415  std::cout<< " The End " << std::endl;
416 
417 
418 }
419 
420 TPZGeoMesh *Hdiv3dPaper201504::GMeshWithPrism( int ndiv)
421 {
422 
423  int Qnodes = 8;
424 
425  TPZGeoMesh * gmesh = new TPZGeoMesh;
426  gmesh->SetMaxNodeId(Qnodes-1);
427  gmesh->NodeVec().Resize(Qnodes);
428  TPZVec<TPZGeoNode> Node(Qnodes);
429 
430  gmesh->SetDimension(3);
431 
432  TPZVec <int64_t> TopolPrism(6);
433  TPZVec <int64_t> TopolQuad(4);
434  TPZVec <int64_t> TopolTriang(3);
435 
436  //indice dos nos
437  int64_t id = 0;
438 
439  TPZManVector<REAL,3> coord(3,0.);
440  int in = 0;
441  //c0
442  coord[0] = 0.0;
443  coord[1] = 0.0;
444  coord[2] = 0.0;
445  gmesh->NodeVec()[in].SetCoord(coord);
446  gmesh->NodeVec()[in].SetNodeId(in);
447  in++;
448  //c1
449  coord[0] = 1.0;
450  coord[1] = 0.0;
451  coord[2] = 0.0;
452  gmesh->NodeVec()[in].SetCoord(coord);
453  gmesh->NodeVec()[in].SetNodeId(in);
454  in++;
455  //c2
456  coord[0] = 1.0;
457  coord[1] = 1.0;
458  coord[2] = 0.0;
459  gmesh->NodeVec()[in].SetCoord(coord);
460  gmesh->NodeVec()[in].SetNodeId(in);
461  in++;
462  //c3
463  coord[0] = 0.0;
464  coord[1] = 1.0;
465  coord[2] = 0.0;
466  gmesh->NodeVec()[in].SetCoord(coord);
467  gmesh->NodeVec()[in].SetNodeId(in);
468  in++;
469  //c4
470  coord[0] = 0.0;
471  coord[1] = 0.0;
472  coord[2] = 1.0;
473  gmesh->NodeVec()[in].SetCoord(coord);
474  gmesh->NodeVec()[in].SetNodeId(in);
475  in++;
476  //c5
477  coord[0] = 1.0;
478  coord[1] = 0.0;
479  coord[2] = 1.0;
480  gmesh->NodeVec()[in].SetCoord(coord);
481  gmesh->NodeVec()[in].SetNodeId(in);
482  in++;
483  //c6
484  coord[0] = 1.0;
485  coord[1] = 1.0;
486  coord[2] = 1.0;
487  gmesh->NodeVec()[in].SetCoord(coord);
488  gmesh->NodeVec()[in].SetNodeId(in);
489  in++;
490  //c7
491  coord[0] = 0.0;
492  coord[1] = 1.0;
493  coord[2] = 1.0;
494  gmesh->NodeVec()[in].SetCoord(coord);
495  gmesh->NodeVec()[in].SetNodeId(in);
496  in++;
497 
498 
499  //indice dos elementos
500  id = 0;
501 
502  TopolTriang[0] = 0;
503  TopolTriang[1] = 1;
504  TopolTriang[2] = 2;
505  new TPZGeoElRefPattern< pzgeom::TPZGeoTriangle > (id,TopolTriang,fbc0,*gmesh);
506  id++;
507 
508  TopolTriang[0] = 0;
509  TopolTriang[1] = 2;
510  TopolTriang[2] = 3;
511  new TPZGeoElRefPattern< pzgeom::TPZGeoTriangle> (id,TopolTriang,fbc0,*gmesh);
512  id++;
513 
514  TopolQuad[0] = 0;
515  TopolQuad[1] = 1;
516  TopolQuad[2] = 5;
517  TopolQuad[3] = 4;
518  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad > (id, TopolQuad, fbc1, *gmesh);
519  id++;
520 
521  TopolQuad[0] = 1;
522  TopolQuad[1] = 2;
523  TopolQuad[2] = 6;
524  TopolQuad[3] = 5;
525  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad > (id, TopolQuad, fbc2, *gmesh);
526  id++;
527 
528  TopolQuad[0] = 3;
529  TopolQuad[1] = 2;
530  TopolQuad[2] = 6;
531  TopolQuad[3] = 7;
532  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad > (id, TopolQuad, fbc3, *gmesh);
533  id++;
534 
535  TopolQuad[0] = 0;
536  TopolQuad[1] = 3;
537  TopolQuad[2] = 7;
538  TopolQuad[3] = 4;
539  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad > (id, TopolQuad, fbc4, *gmesh);
540  id++;
541 
542  TopolTriang[0] = 4;
543  TopolTriang[1] = 5;
544  TopolTriang[2] = 6;
545  new TPZGeoElRefPattern< pzgeom::TPZGeoTriangle > (id,TopolTriang,fbc5,*gmesh);
546  id++;
547 
548  TopolTriang[0] = 4;
549  TopolTriang[1] = 6;
550  TopolTriang[2] = 7;
551  new TPZGeoElRefPattern< pzgeom::TPZGeoTriangle> (id,TopolTriang,fbc5,*gmesh);
552  id++;
553 
554  TopolPrism[0] = 0;
555  TopolPrism[1] = 1;
556  TopolPrism[2] = 2;
557  TopolPrism[3] = 4;
558  TopolPrism[4] = 5;
559  TopolPrism[5] = 6;
560  new TPZGeoElRefPattern< pzgeom::TPZGeoPrism > (id, TopolPrism,fmatId,*gmesh);
561  id++;
562 
563  TopolPrism[0] = 0;
564  TopolPrism[1] = 2;
565  TopolPrism[2] = 3;
566  TopolPrism[3] = 4;
567  TopolPrism[4] = 6;
568  TopolPrism[5] = 7;
569  new TPZGeoElRefPattern< pzgeom::TPZGeoPrism > (id, TopolPrism,fmatId,*gmesh);
570  id++;
571 
572 
573  gmesh->BuildConnectivity();
574 
576 
577  TPZVec<TPZGeoEl *> sons;
578  for (int iref = 0; iref < ndiv; iref++) {
579  int nel = gmesh->NElements();
580  for (int iel = 0; iel < nel; iel++) {
581  TPZGeoEl *gel = gmesh->ElementVec()[iel];
582  if (gel->HasSubElement()) {
583  continue;
584  }
585  gel->Divide(sons);
586  }
587  }
588 
589 
590  //#ifdef LOG4CXX
591  // if(logdata->isDebugEnabled())
592  // {
593  // std::stringstream sout;
594  // sout<<"\n\n Malha Geometrica Inicial\n ";
595  // gmesh->Print(sout);
596  // LOGPZ_DEBUG(logdata,sout.str())
597  // }
598  //#endif
599 
600  return gmesh;
601 }
602 
603 TPZGeoMesh *Hdiv3dPaper201504::CreateOneCuboWithTetraedrons(int64_t nelem)
604 {
605  TPZGeoMesh *gmesh = new TPZGeoMesh;
606  GenerateNodes(gmesh,nelem);
607 
608  for (int64_t i=0; i<nelem; i++) {
609  for (int64_t j=0; j<nelem; j++) {
610  for (int64_t k=0; k<nelem; k++) {
611  TPZManVector<int64_t,8> nodes(8,0);
612  nodes[0] = k*(nelem+1)*(nelem+1)+j*(nelem+1)+i;
613  nodes[1] = k*(nelem+1)*(nelem+1)+j*(nelem+1)+i+1;
614  nodes[2] = k*(nelem+1)*(nelem+1)+(j+1)*(nelem+1)+i+1;
615  nodes[3] = k*(nelem+1)*(nelem+1)+(j+1)*(nelem+1)+i;
616  nodes[4] = (k+1)*(nelem+1)*(nelem+1)+j*(nelem+1)+i;
617  nodes[5] = (k+1)*(nelem+1)*(nelem+1)+j*(nelem+1)+i+1;
618  nodes[6] = (k+1)*(nelem+1)*(nelem+1)+(j+1)*(nelem+1)+i+1;
619  nodes[7] = (k+1)*(nelem+1)*(nelem+1)+(j+1)*(nelem+1)+i;
620 
621  for (int el=0; el<6; el++)
622  {
623  TPZManVector<int64_t,4> elnodes(4);
624  int64_t index;
625  for (int il=0; il<4; il++) {
626  elnodes[il] = nodes[tetraedra_2(el,il)];//nodes[tetraedra_2[el][il]];
627  }
628  gmesh->CreateGeoElement(ETetraedro, elnodes, fmatId, index);
629  }
630  }
631  }
632  }
633 
634  gmesh->BuildConnectivity();
635 
636  // Boundary Conditions
637  const int numelements = gmesh->NElements();
638  // const int bczMinus = -3, bczplus = -2, bcids = -1;
639  // const int bczMinus = -1, bczplus = -1, bcids = -1;
640 
641  for(int el=0; el<numelements; el++)
642  {
643  TPZManVector <TPZGeoNode,4> Nodefinder(4);
644  TPZManVector <REAL,3> nodecoord(3);
645  TPZGeoEl *tetra = gmesh->ElementVec()[el];
646  TPZVec<int64_t> ncoordVec(0);
647  int64_t sizeOfVec = 0;
648 
649  // na face z = 0
650  for (int i = 0; i < 4; i++)
651  {
652  int64_t pos = tetra->NodeIndex(i);
653  Nodefinder[i] = gmesh->NodeVec()[pos];
654  Nodefinder[i].GetCoordinates(nodecoord);
655  if (MyDoubleComparer(nodecoord[2],0.))
656  {
657  sizeOfVec++;
658  ncoordVec.Resize(sizeOfVec);
659  ncoordVec[sizeOfVec-1] = pos;
660  }
661  }
662  if(ncoordVec.NElements() == 3)
663  {
664  int lado = tetra->WhichSide(ncoordVec);
665  TPZGeoElSide tetraSide(tetra, lado);
666  TPZGeoElBC(tetraSide,fbc0);
667  }
668 
669  ncoordVec.clear();
670  sizeOfVec = 0;
671  // na face y = 0
672  for (int i = 0; i < 4; i++)
673  {
674  int64_t pos = tetra->NodeIndex(i);
675  Nodefinder[i] = gmesh->NodeVec()[pos];
676  Nodefinder[i].GetCoordinates(nodecoord);
677  if (MyDoubleComparer(nodecoord[1],0.))
678  {
679  sizeOfVec++;
680  ncoordVec.Resize(sizeOfVec);
681  ncoordVec[sizeOfVec-1] = pos;
682  }
683  }
684  if(ncoordVec.NElements() == 3)
685  {
686  int lado = tetra->WhichSide(ncoordVec);
687  TPZGeoElSide tetraSide(tetra, lado);
688  TPZGeoElBC(tetraSide,fbc1);
689  }
690 
691  ncoordVec.clear();
692  sizeOfVec = 0;
693  // na face x = 1
694  for (int i = 0; i < 4; i++)
695  {
696  int64_t pos = tetra->NodeIndex(i);
697  Nodefinder[i] = gmesh->NodeVec()[pos];
698  Nodefinder[i].GetCoordinates(nodecoord);
699  if (MyDoubleComparer(nodecoord[0],1.))
700  {
701  sizeOfVec++;
702  ncoordVec.Resize(sizeOfVec);
703  ncoordVec[sizeOfVec-1] = pos;
704  }
705  }
706  if(ncoordVec.NElements() == 3)
707  {
708  int lado = tetra->WhichSide(ncoordVec);
709  TPZGeoElSide tetraSide(tetra, lado);
710  TPZGeoElBC(tetraSide,fbc2);
711  }
712 
713  ncoordVec.clear();
714  sizeOfVec = 0;
715  // na face y = 1
716  for (int i = 0; i < 4; i++)
717  {
718  int64_t pos = tetra->NodeIndex(i);
719  Nodefinder[i] = gmesh->NodeVec()[pos];
720  Nodefinder[i].GetCoordinates(nodecoord);
721  if (MyDoubleComparer(nodecoord[1],1.))
722  {
723  sizeOfVec++;
724  ncoordVec.Resize(sizeOfVec);
725  ncoordVec[sizeOfVec-1] = pos;
726  }
727  }
728  if(ncoordVec.NElements() == 3)
729  {
730  int lado = tetra->WhichSide(ncoordVec);
731  TPZGeoElSide tetraSide(tetra, lado);
732  TPZGeoElBC(tetraSide,fbc3);
733  }
734 
735 
736  ncoordVec.clear();
737  sizeOfVec = 0;
738  // na face x = 0
739  for (int i = 0; i < 4; i++)
740  {
741  int64_t pos = tetra->NodeIndex(i);
742  Nodefinder[i] = gmesh->NodeVec()[pos];
743  Nodefinder[i].GetCoordinates(nodecoord);
744  if (MyDoubleComparer(nodecoord[0],0.))
745  {
746  sizeOfVec++;
747  ncoordVec.Resize(sizeOfVec);
748  ncoordVec[sizeOfVec-1] = pos;
749  }
750  }
751  if(ncoordVec.NElements() == 3)
752  {
753  int lado = tetra->WhichSide(ncoordVec);
754  TPZGeoElSide tetraSide(tetra, lado);
755  TPZGeoElBC(tetraSide,fbc4);
756  }
757 
758  ncoordVec.clear();
759  sizeOfVec = 0;
760  // na face z = 1
761  for (int i = 0; i < 4; i++)
762  {
763  int64_t pos = tetra->NodeIndex(i);
764  Nodefinder[i] = gmesh->NodeVec()[pos];
765  Nodefinder[i].GetCoordinates(nodecoord);
766  if (MyDoubleComparer(nodecoord[2],1.))
767  {
768  sizeOfVec++;
769  ncoordVec.Resize(sizeOfVec);
770  ncoordVec[sizeOfVec-1] = pos;
771  }
772  }
773  if(ncoordVec.NElements() == 3)
774  {
775  int lado = tetra->WhichSide(ncoordVec);
776  TPZGeoElSide tetraSide(tetra, lado);
777  TPZGeoElBC(tetraSide,fbc5);
778  }
779 
780 
781 
782  }
783 
784  return gmesh;
785 }
786 
787 void Hdiv3dPaper201504::GenerateNodes(TPZGeoMesh *gmesh, int64_t nelem)
788 {
789  gmesh->NodeVec().Resize((nelem+1)*(nelem+1)*(nelem+1));
790  for (int64_t i=0; i<=nelem; i++) {
791  for (int64_t j=0; j<=nelem; j++) {
792  for (int64_t k=0; k<=nelem; k++) {
793  TPZManVector<REAL,3> x(3);
794  x[0] = k*1./nelem;
795  x[1] = j*1./nelem;
796  x[2] = i*1./nelem;
797  gmesh->NodeVec()[i*(nelem+1)*(nelem+1)+j*(nelem+1)+k].Initialize(x, *gmesh);
798  }
799  }
800  }
801 }
802 
803 
804 TPZGeoMesh *Hdiv3dPaper201504::CreateOneCubo(int nref)
805 {
806  TPZGeoMesh *gmesh = new TPZGeoMesh;
807  int nnodes = 8;
808 
809  gmesh->SetDimension(3);
810  gmesh->NodeVec().Resize(nnodes);
811 
812  TPZManVector<REAL,3> coord(3,0.);
813  int in = 0;
814 
815  //cubo [0,1]ˆ3
816  //c0
817  coord[0] = 0.0;
818  coord[1] = 0.0;
819  coord[2] = 0.0;
820  gmesh->NodeVec()[in].SetCoord(coord);
821  gmesh->NodeVec()[in].SetNodeId(in);
822  in++;
823  //c1
824  coord[0] = 1.0;
825  coord[1] = 0.0;
826  coord[2] = 0.0;
827  gmesh->NodeVec()[in].SetCoord(coord);
828  gmesh->NodeVec()[in].SetNodeId(in);
829  in++;
830  //c2
831  coord[0] = 1.0;
832  coord[1] = 1.0;
833  coord[2] = 0.0;
834  gmesh->NodeVec()[in].SetCoord(coord);
835  gmesh->NodeVec()[in].SetNodeId(in);
836  in++;
837  //c3
838  coord[0] = 0.0;
839  coord[1] = 1.0;
840  coord[2] = 0.0;
841  gmesh->NodeVec()[in].SetCoord(coord);
842  gmesh->NodeVec()[in].SetNodeId(in);
843  in++;
844 
845  //c4
846  coord[0] = 0.0;
847  coord[1] = 0.0;
848  coord[2] = 1.0;
849  gmesh->NodeVec()[in].SetCoord(coord);
850  gmesh->NodeVec()[in].SetNodeId(in);
851  in++;
852  //c5
853  coord[0] = 1.0;
854  coord[1] = 0.0;
855  coord[2] = 1.0;
856  gmesh->NodeVec()[in].SetCoord(coord);
857  gmesh->NodeVec()[in].SetNodeId(in);
858  in++;
859  //c6
860  coord[0] = 1.0;
861  coord[1] = 1.0;
862  coord[2] = 1.0;
863  gmesh->NodeVec()[in].SetCoord(coord);
864  gmesh->NodeVec()[in].SetNodeId(in);
865  in++;
866  //c7
867  coord[0] = 0.0;
868  coord[1] = 1.0;
869  coord[2] = 1.0;
870  gmesh->NodeVec()[in].SetCoord(coord);
871  gmesh->NodeVec()[in].SetNodeId(in);
872  in++;
873 
874  // cubo [-1,1]^3
875  // //c0
876  // coord[0] = -1.0;
877  // coord[1] = -1.0;
878  // coord[2] = -1.0;
879  // gmesh->NodeVec()[in].SetCoord(coord);
880  // gmesh->NodeVec()[in].SetNodeId(in);
881  // in++;
882  // //c1
883  // coord[0] = 1.0;
884  // coord[1] = -1.0;
885  // coord[2] = -1.0;
886  // gmesh->NodeVec()[in].SetCoord(coord);
887  // gmesh->NodeVec()[in].SetNodeId(in);
888  // in++;
889  // //c2
890  // coord[0] = 1.0;
891  // coord[1] = 1.0;
892  // coord[2] = -1.0;
893  // gmesh->NodeVec()[in].SetCoord(coord);
894  // gmesh->NodeVec()[in].SetNodeId(in);
895  // in++;
896  // //c3
897  // coord[0] = -1.0;
898  // coord[1] = 1.0;
899  // coord[2] = -1.0;
900  // gmesh->NodeVec()[in].SetCoord(coord);
901  // gmesh->NodeVec()[in].SetNodeId(in);
902  // in++;
903  // //c4
904  // coord[0] = -1.0;
905  // coord[1] = -1.0;
906  // coord[2] = 1.0;
907  // gmesh->NodeVec()[in].SetCoord(coord);
908  // gmesh->NodeVec()[in].SetNodeId(in);
909  // in++;
910  // //c5
911  // coord[0] = 1.0;
912  // coord[1] = -1.0;
913  // coord[2] = 1.0;
914  // gmesh->NodeVec()[in].SetCoord(coord);
915  // gmesh->NodeVec()[in].SetNodeId(in);
916  // in++;
917  // //c6
918  // coord[0] = 1.0;
919  // coord[1] = 1.0;
920  // coord[2] = 1.0;
921  // gmesh->NodeVec()[in].SetCoord(coord);
922  // gmesh->NodeVec()[in].SetNodeId(in);
923  // in++;
924  // //c7
925  // coord[0] = -1.0;
926  // coord[1] = 1.0;
927  // coord[2] = 1.0;
928  // gmesh->NodeVec()[in].SetCoord(coord);
929  // gmesh->NodeVec()[in].SetNodeId(in);
930  // in++;
931 
932 
933 
934  int index = 0;
935 
936  TPZVec<int64_t> TopologyQuad(4);
937 
938  // bottom
939  TopologyQuad[0] = 0;
940  TopologyQuad[1] = 1;
941  TopologyQuad[2] = 2;
942  TopologyQuad[3] = 3;
943  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc0,*gmesh);
944  index++;
945 
946  // Front
947  TopologyQuad[0] = 0;
948  TopologyQuad[1] = 1;
949  TopologyQuad[2] = 5;
950  TopologyQuad[3] = 4;
951  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc1,*gmesh);
952  index++;
953 
954  // Rigth
955  TopologyQuad[0] = 1;
956  TopologyQuad[1] = 2;
957  TopologyQuad[2] = 6;
958  TopologyQuad[3] = 5;
959  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc2,*gmesh);
960  index++;
961  // Back
962  TopologyQuad[0] = 3;
963  TopologyQuad[1] = 2;
964  TopologyQuad[2] = 6;
965  TopologyQuad[3] = 7;
966  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc3,*gmesh);
967  index++;
968 
969  // Left
970  TopologyQuad[0] = 0;
971  TopologyQuad[1] = 3;
972  TopologyQuad[2] = 7;
973  TopologyQuad[3] = 4;
974  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc4,*gmesh);
975  index++;
976 
977  // Top
978  TopologyQuad[0] = 4;
979  TopologyQuad[1] = 5;
980  TopologyQuad[2] = 6;
981  TopologyQuad[3] = 7;
982  new TPZGeoElRefPattern< pzgeom::TPZGeoQuad>(index,TopologyQuad,fbc5,*gmesh);
983  index++;
984 
985  TPZManVector<int64_t,8> TopolCubo(8,0);
986  TopolCubo[0] = 0;
987  TopolCubo[1] = 1;
988  TopolCubo[2] = 2;
989  TopolCubo[3] = 3;
990  TopolCubo[4] = 4;
991  TopolCubo[5] = 5;
992  TopolCubo[6] = 6;
993  TopolCubo[7] = 7;
994 
995 
996  new TPZGeoElRefPattern< pzgeom::TPZGeoCube> (index, TopolCubo, fmatId, *gmesh);
997  index++;
998 
999 
1000  gmesh->BuildConnectivity();
1001 
1003 
1004  TPZVec<TPZGeoEl *> sons;
1005  for (int iref = 0; iref < nref; iref++) {
1006  int nel = gmesh->NElements();
1007  for (int iel = 0; iel < nel; iel++) {
1008  TPZGeoEl *gel = gmesh->ElementVec()[iel];
1009  if (gel->HasSubElement()) {
1010  continue;
1011  }
1012  gel->Divide(sons);
1013  }
1014  }
1015 
1016 
1017  std::ofstream out("SingleCubeWithBcs.vtk");
1018  TPZVTKGeoMesh::PrintGMeshVTK(gmesh, out, true);
1019 
1020  return gmesh;
1021 }
1022 
1023 
1024 
1025 void Hdiv3dPaper201504::SolExata(const TPZVec<REAL> &pt, TPZVec<STATE> &solp, TPZFMatrix<STATE> &flux){
1026 
1027  solp.resize(1);
1028  solp[0]=0.;
1029  flux.Resize(3, 1);
1030  double x = pt[0];
1031  double y = pt[1];
1032  double z = pt[2];
1033  for(int d=0; d<3;d++)
1034  {
1035  flux(d,0)=0.;
1036  }
1037 
1038 
1039  solp[0] = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1040  flux(0,0) = -M_PI*(cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z));
1041  flux(1,0) = -M_PI*(sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z));
1042  flux(2,0) = -M_PI*(sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z));
1043 
1044 }
1045 
1046 void Hdiv3dPaper201504::Forcing(const TPZVec<REAL> &pt, TPZVec<STATE> &ff){
1047 
1048  double x = pt[0];
1049  double y = pt[1];
1050  double z = pt[2];
1051 
1052  ff[0] = 3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1053 
1054 }
1055 
1056 void Hdiv3dPaper201504::SolExataH1(const TPZVec<REAL> &pt, TPZVec<STATE> &solp, TPZFMatrix<STATE> &flux){
1057 
1058  solp.resize(1);
1059  solp[0]=0.;
1060  flux.Resize(3, 1);
1061  double x = pt[0];
1062  double y = pt[1];
1063  double z = pt[2];
1064  for(int d=0; d<3;d++)
1065  {
1066  flux(d,0)=0.;
1067  }
1068 
1069 
1070  solp[0] = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1071  flux(0,0) = M_PI*(cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z));
1072  flux(1,0) = M_PI*(sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z));
1073  flux(2,0) = M_PI*(sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z));
1074 
1075 }
1076 
1077 void Hdiv3dPaper201504::ForcingH1(const TPZVec<REAL> &pt, TPZVec<STATE> &ff, TPZFMatrix<STATE> &flux)
1078 {
1079 
1080  double x = pt[0];
1081  double y = pt[1];
1082  double z = pt[2];
1083 
1084  ff[0] = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1085 
1086 }
1087 
1088 
1089 void Hdiv3dPaper201504::ForcingBC0D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1090 
1091  solp[0] = 0.0;
1092 
1093 }
1094 
1095 void Hdiv3dPaper201504::ForcingBC1D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1096 
1097  solp[0] = 0.0;
1098 
1099 }
1100 
1101 void Hdiv3dPaper201504::ForcingBC2D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1102 
1103  solp[0] = 0.0;
1104 
1105 }
1106 
1107 void Hdiv3dPaper201504::ForcingBC3D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1108 
1109  solp[0] = 0.0;
1110 
1111 }
1112 
1113 void Hdiv3dPaper201504::ForcingBC4D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1114 
1115 
1116  solp[0] = 0.0;
1117 
1118 }
1119 
1120 void Hdiv3dPaper201504::ForcingBC5D(const TPZVec<REAL> &pt, TPZVec<STATE> &solp){
1121 
1122  solp[0] = 0.0;
1123 
1124 }
1125 
1126 TPZCompMesh *Hdiv3dPaper201504::CMeshH1(TPZGeoMesh *gmesh, int pOrder, int dim)
1127 {
1129  TPZMatPoisson3d *material = new TPZMatPoisson3d(fmatId,fDim);
1130  TPZMaterial * mat(material);
1131  material->NStateVariables();
1132 
1133  // //solucao exata
1134  TPZAutoPointer<TPZFunction<STATE> > solexata;
1135  solexata = new TPZDummyFunction<STATE>(SolExataH1, 5);
1136  material->SetForcingFunctionExact(solexata);
1137 
1138  //funcao do lado direito da equacao do problema
1139  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(ForcingH1, 5);
1140  TPZAutoPointer<TPZFunction<STATE> > forcef;
1141  dum->SetPolynomialOrder(20);
1142  forcef = dum;
1143  material->SetForcingFunction(forcef);
1144 
1145  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1146  cmesh->SetDimModel(fDim);
1147  cmesh->SetDefaultOrder(pOrder);
1148 
1149  cmesh->InsertMaterialObject(mat);
1150 
1152  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1153  val2(0,0) = 0.0;
1154  val2(1,0) = 0.0;
1155  TPZMaterial * BCond0;
1156  TPZMaterial * BCond1;
1157  TPZMaterial * BCond2;
1158  TPZMaterial * BCond3;
1159  TPZMaterial * BCond4;
1160  TPZMaterial * BCond5;
1161 
1162  BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);
1163  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1164  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1165  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1166  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1167  BCond5 = material->CreateBC(mat, fbc5,fdirichlet, val1, val2);
1168 
1169  cmesh->InsertMaterialObject(BCond0);
1170  cmesh->InsertMaterialObject(BCond1);
1171  cmesh->InsertMaterialObject(BCond2);
1172  cmesh->InsertMaterialObject(BCond3);
1173  cmesh->InsertMaterialObject(BCond4);
1174  cmesh->InsertMaterialObject(BCond5);
1175 
1176 
1177  cmesh->SetAllCreateFunctionsContinuous();
1178 
1179  //Ajuste da estrutura de dados computacional
1180  cmesh->AutoBuild();
1181 
1182  return cmesh;
1183 
1184 }
1185 
1186 
1187 TPZCompMesh *Hdiv3dPaper201504::CMeshFlux(TPZGeoMesh *gmesh, int pOrder, int dim)
1188 {
1190  //TPZMatPoisson3d *material = new TPZMatPoisson3d(matId,fDim);
1191  //TPZMatPoissonD3 *material = new TPZMatPoissonD3(fmatId,fDim);
1192  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
1193 
1194  TPZMaterial * mat(material);
1195  material->NStateVariables();
1196 
1197  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1198 
1199 
1201  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1202  TPZMaterial * BCond0;
1203  TPZMaterial * BCond1;
1204  TPZMaterial * BCond2;
1205  TPZMaterial * BCond3;
1206  TPZMaterial * BCond4;
1207  TPZMaterial * BCond5;
1208 
1209  BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);
1210  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1211  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1212  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1213  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1214  BCond5 = material->CreateBC(mat, fbc5,fdirichlet, val1, val2);
1215 
1216  cmesh->InsertMaterialObject(mat);
1217 
1218  cmesh->SetDimModel(fDim);
1219 
1220  cmesh->SetAllCreateFunctionsHDiv();
1221 
1222  cmesh->InsertMaterialObject(BCond0);
1223  cmesh->InsertMaterialObject(BCond1);
1224  cmesh->InsertMaterialObject(BCond2);
1225  cmesh->InsertMaterialObject(BCond3);
1226  cmesh->InsertMaterialObject(BCond4);
1227  cmesh->InsertMaterialObject(BCond5);
1228 
1229 
1230  cmesh->SetDefaultOrder(pOrder);
1231 
1232 
1233  TPZLagrangeMultiplier *matskelet = new TPZLagrangeMultiplier(fmatskeleton, fDim-1, 1);
1234  TPZMaterial * mat2(matskelet);
1235  cmesh->InsertMaterialObject(mat2);
1236 
1237 
1238 
1239  //Ajuste da estrutura de dados computacional
1240  cmesh->AutoBuild();
1241 
1242  return cmesh;
1243 
1244 }
1245 
1246 TPZCompMesh *Hdiv3dPaper201504::CMeshPressure(TPZGeoMesh *gmesh, int pOrder, int dim)
1247 {
1249  //TPZMatPoisson3d *material = new TPZMatPoisson3d(fmatId,fDim);
1250  //TPZMatPoissonD3 *material = new TPZMatPoissonD3(fmatId,fDim);
1251  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,fDim);
1252  material->NStateVariables();
1253 
1254  TPZCompMesh * cmesh = new TPZCompMesh(gmesh);
1255  cmesh->SetDimModel(fDim);
1256  TPZMaterial * mat(material);
1257  cmesh->InsertMaterialObject(mat);
1258 
1260  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1261 
1262  cmesh->SetDefaultOrder(pOrder);
1263 
1264  bool h1function = true;// com esqueleto precisa disso
1265  if(pOrder>0/*h1function*/){
1266  cmesh->SetAllCreateFunctionsContinuous();
1267  cmesh->ApproxSpace().CreateDisconnectedElements(true);
1268  }
1269  else{
1270  cmesh->SetAllCreateFunctionsDiscontinuous();
1271  }
1272 
1273 
1274  //Ajuste da estrutura de dados computacional
1275  cmesh->AutoBuild();
1276 
1277 
1278  int ncon = cmesh->NConnects();
1279  for(int i=0; i<ncon; i++)
1280  {
1281  TPZConnect &newnod = cmesh->ConnectVec()[i];
1282  //newnod.SetPressure(true);
1283  newnod.SetLagrangeMultiplier(1);
1284  }
1285 
1286  if(!h1function)
1287  {
1288 
1289  int nel = cmesh->NElements();
1290  for(int i=0; i<nel; i++){
1291  TPZCompEl *cel = cmesh->ElementVec()[i];
1292  TPZCompElDisc *celdisc = dynamic_cast<TPZCompElDisc *>(cel);
1293  celdisc->SetConstC(1.);
1294  celdisc->SetCenterPoint(0, 0.);
1295  celdisc->SetCenterPoint(1, 0.);
1296  celdisc->SetCenterPoint(2, 0.);
1297  celdisc->SetTrueUseQsiEta();
1298 
1299  if(celdisc && celdisc->Reference()->Dimension() == cmesh->Dimension())
1300  {
1301  //if(ftriang==true) celdisc->SetTotalOrderShape();
1302  //else celdisc->SetTensorialShape();
1303  }
1304 
1305  }
1306  }
1307 
1308 #ifdef PZDEBUG
1309  int ncel = cmesh->NElements();
1310  for(int i =0; i<ncel; i++){
1311  TPZCompEl * compEl = cmesh->ElementVec()[i];
1312  if(!compEl) continue;
1313  TPZInterfaceElement * facel = dynamic_cast<TPZInterfaceElement *>(compEl);
1314  if(facel)DebugStop();
1315 
1316  }
1317 #endif
1318 
1319  return cmesh;
1320 
1321 }
1322 
1323 TPZCompMesh *Hdiv3dPaper201504::CMeshMixed(TPZGeoMesh * gmesh, TPZVec<TPZCompMesh *> meshvec)
1324 {
1325 
1326  //Creating computational mesh for multiphysic elements
1327  gmesh->ResetReference();
1328  TPZCompMesh *mphysics = new TPZCompMesh(gmesh);
1329 
1330  //criando material
1331  int dim = gmesh->Dimension();
1332 
1333  TPZMatMixedPoisson3D *material = new TPZMatMixedPoisson3D(fmatId,dim);
1334 
1335  //incluindo os dados do problema
1336  TPZFNMatrix<9,REAL> PermTensor(3,3,0.);
1337  TPZFNMatrix<9,REAL> InvPermTensor(3,3,0.);
1338 
1339 
1340  // tensor de permutacao
1341  TPZFNMatrix<9,REAL> TP(3,3,0.0);
1342  TPZFNMatrix<9,REAL> InvTP(3,3,0.0);
1343 
1344  // Hard coded
1345  for (int id = 0; id < 3; id++){
1346  TP(id,id) = 1.0;
1347  InvTP(id,id) = 1.0;
1348  }
1349 
1350  PermTensor = TP;
1351  InvPermTensor = InvTP;
1352 
1353  material->SetPermeabilityTensor(PermTensor, InvPermTensor);
1354 
1355  //solucao exata
1356  TPZAutoPointer<TPZFunction<STATE> > solexata;
1357 
1358  solexata = new TPZDummyFunction<STATE>(SolExata, 5);
1359  material->SetForcingFunctionExact(solexata);
1360  mphysics->SetDimModel(dim);
1361  //funcao do lado direito da equacao do problema
1362  TPZDummyFunction<STATE> *dum = new TPZDummyFunction<STATE>(Forcing, 5);
1363  TPZAutoPointer<TPZFunction<STATE> > forcef;
1364  dum->SetPolynomialOrder(10);
1365  forcef = dum;
1366  material->SetForcingFunction(forcef);
1367 
1368  //inserindo o material na malha computacional
1369  TPZMaterial *mat(material);
1370  mphysics->InsertMaterialObject(mat);
1371 
1372 
1373  //Criando condicoes de contorno
1374  TPZMaterial * BCond0;
1375  TPZMaterial * BCond1;
1376  TPZMaterial * BCond2;
1377  TPZMaterial * BCond3;
1378  TPZMaterial * BCond4;
1379  TPZMaterial * BCond5;
1380 
1381  TPZFMatrix<STATE> val1(2,2,0.), val2(2,1,0.);
1382 
1383  val2(0,0) = 0.0;
1384  val2(1,0) = 0.0;
1385  TPZAutoPointer<TPZFunction<STATE> > FBCond0 = new TPZDummyFunction<STATE>(ForcingBC0D, 5);
1386  BCond0 = material->CreateBC(mat, fbc0,fdirichlet, val1, val2);
1387  BCond0->SetForcingFunction(FBCond0);
1388 
1389  val2(0,0) = 0.0;
1390  val2(1,0) = 0.0;
1391  TPZAutoPointer<TPZFunction<STATE> > FBCond1 = new TPZDummyFunction<STATE>(ForcingBC1D, 5);
1392  BCond1 = material->CreateBC(mat, fbc1,fdirichlet, val1, val2);
1393  BCond1->SetForcingFunction(FBCond1);
1394 
1395  val2(0,0) = 0.0;
1396  val2(1,0) = 0.0;
1397  TPZAutoPointer<TPZFunction<STATE> > FBCond2 = new TPZDummyFunction<STATE>(ForcingBC2D, 5);
1398  BCond2 = material->CreateBC(mat, fbc2,fdirichlet, val1, val2);
1399  BCond2->SetForcingFunction(FBCond2);
1400 
1401  val2(0,0) = 0.0;
1402  val2(1,0) = 0.0;
1403  TPZAutoPointer<TPZFunction<STATE> > FBCond3 = new TPZDummyFunction<STATE>(ForcingBC3D, 5);
1404  BCond3 = material->CreateBC(mat, fbc3,fdirichlet, val1, val2);
1405  BCond3->SetForcingFunction(FBCond3);
1406 
1407  val2(0,0) = 0.0;
1408  val2(1,0) = 0.0;
1409  TPZAutoPointer<TPZFunction<STATE> > FBCond4 = new TPZDummyFunction<STATE>(ForcingBC4D, 5);
1410  BCond4 = material->CreateBC(mat, fbc4,fdirichlet, val1, val2);
1411  BCond4->SetForcingFunction(FBCond4);
1412 
1413  val2(0,0) = 0.0;
1414  val2(1,0) = 0.0;
1415  TPZAutoPointer<TPZFunction<STATE> > FBCond5 = new TPZDummyFunction<STATE>(ForcingBC5D, 5);
1416  BCond5 = material->CreateBC(mat, fbc5,fdirichlet, val1, val2);
1417  BCond5->SetForcingFunction(FBCond5);
1418 
1419 
1420 
1421  mphysics->SetAllCreateFunctionsMultiphysicElem();
1422 
1423  mphysics->InsertMaterialObject(BCond0);
1424  mphysics->InsertMaterialObject(BCond1);
1425  mphysics->InsertMaterialObject(BCond2);
1426  mphysics->InsertMaterialObject(BCond3);
1427  mphysics->InsertMaterialObject(BCond4);
1428  mphysics->InsertMaterialObject(BCond5);
1429 
1430  //Fazendo auto build
1431  mphysics->AutoBuild();
1432  mphysics->AdjustBoundaryElements();
1433  mphysics->CleanUpUnconnectedNodes();
1434 
1435  //Creating multiphysic elements containing skeletal elements.
1436  TPZBuildMultiphysicsMesh::AddElements(meshvec, mphysics);
1437  mphysics->Reference()->ResetReference();
1438  mphysics->LoadReferences();
1439 
1440  int64_t nel = mphysics->ElementVec().NElements();
1441 
1442  std::map<int64_t, int64_t> bctoel, eltowrap;
1443  for (int64_t el=0; el<nel; el++) {
1444  TPZCompEl *cel = mphysics->Element(el);
1445  TPZGeoEl *gel = cel->Reference();
1446  int matid = gel->MaterialId();
1447  if (matid < 0) {
1448  TPZGeoElSide gelside(gel,gel->NSides()-1);
1449  TPZGeoElSide neighbour = gelside.Neighbour();
1450  while (neighbour != gelside) {
1451  if (neighbour.Element()->Dimension() == dim && neighbour.Element()->Reference()) {
1452  // got you!!
1453  bctoel[el] = neighbour.Element()->Reference()->Index();
1454  break;
1455  }
1456  neighbour = neighbour.Neighbour();
1457  }
1458  if (neighbour == gelside) {
1459  DebugStop();
1460  }
1461  }
1462  }
1463 
1464  TPZStack< TPZStack< TPZMultiphysicsElement *,7> > wrapEl;
1465  for(int64_t el = 0; el < nel; el++)
1466  {
1467  TPZMultiphysicsElement *mfcel = dynamic_cast<TPZMultiphysicsElement *>(mphysics->Element(el));
1468  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
1469  }
1470 
1471  for (int64_t el =0; el < wrapEl.size(); el++) {
1472  TPZCompEl *cel = wrapEl[el][0];
1473  int64_t index = cel->Index();
1474  eltowrap[index] = el;
1475  }
1476 
1477  meshvec[0]->CleanUpUnconnectedNodes();
1478  TPZBuildMultiphysicsMesh::AddConnects(meshvec,mphysics);
1479  TPZBuildMultiphysicsMesh::TransferFromMeshes(meshvec, mphysics);
1480 
1481  std::map<int64_t, int64_t>::iterator it;
1482  for (it = bctoel.begin(); it != bctoel.end(); it++) {
1483  int64_t bcindex = it->first;
1484  int64_t elindex = it->second;
1485  if (eltowrap.find(elindex) == eltowrap.end()) {
1486  DebugStop();
1487  }
1488  int64_t wrapindex = eltowrap[elindex];
1489  TPZCompEl *bcel = mphysics->Element(bcindex);
1490  TPZMultiphysicsElement *bcmf = dynamic_cast<TPZMultiphysicsElement *>(bcel);
1491  if (!bcmf) {
1492  DebugStop();
1493  }
1494  wrapEl[wrapindex].Push(bcmf);
1495 
1496  }
1497 
1498  //------- Create and add group elements -------
1499  int64_t index, nenvel;
1500  nenvel = wrapEl.NElements();
1501  TPZStack<TPZElementGroup *> elgroups;
1502  for(int64_t ienv=0; ienv<nenvel; ienv++){
1503  TPZElementGroup *elgr = new TPZElementGroup(*wrapEl[ienv][0]->Mesh(),index);
1504  elgroups.Push(elgr);
1505  nel = wrapEl[ienv].NElements();
1506  for(int jel=0; jel<nel; jel++){
1507  elgr->AddElement(wrapEl[ienv][jel]);
1508  }
1509  }
1510 
1511  mphysics->ComputeNodElCon();
1512  // create condensed elements
1513  // increase the NumElConnected of one pressure connects in order to prevent condensation
1514  for (int64_t ienv=0; ienv<nenvel; ienv++) {
1515  TPZElementGroup *elgr = elgroups[ienv];
1516  int nc = elgr->NConnects();
1517  for (int ic=0; ic<nc; ic++) {
1518  TPZConnect &c = elgr->Connect(ic);
1519  if (c.LagrangeMultiplier() > 0) {
1520  c.IncrementElConnected();
1521  break;
1522  }
1523  }
1524  new TPZCondensedCompEl(elgr);
1525  }
1526 
1527  mphysics->CleanUpUnconnectedNodes();
1528  mphysics->ExpandSolution();
1529 
1530  return mphysics;
1531 
1532 }
1533 
1534 
1535 void Hdiv3dPaper201504::ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
1536 {
1537 
1538  int64_t nel = l2mesh->NElements();
1539  int dim = l2mesh->Dimension();
1540  TPZManVector<STATE,10> globalerrors(10,0.);
1541  for (int64_t el=0; el<nel; el++) {
1542  TPZCompEl *cel = l2mesh->ElementVec()[el];
1543  if (!cel) {
1544  continue;
1545  }
1546  TPZGeoEl *gel = cel->Reference();
1547  if (!gel || gel->Dimension() != dim) {
1548  continue;
1549  }
1550  TPZManVector<REAL,10> elerror(10,0.);
1551  elerror.Fill(0.);
1552  cel->EvaluateError(SolExataH1, elerror, 0);
1553 
1554  int nerr = elerror.size();
1555  globalerrors.resize(nerr);
1556 
1557  for (int i=0; i<nerr; i++) {
1558  globalerrors[i] += elerror[i]*elerror[i];
1559  }
1560  }
1561  // sequence of storage
1562  // for each
1563  errors.PutVal(pos, 0, p); // polinomial order
1564  errors.PutVal(pos, 1, ndiv); // refinement
1565  errors.PutVal(pos, 2, sqrt(globalerrors[1])); // pressure
1566  errors.PutVal(pos, 3, sqrt(globalerrors[2])); // flux
1567 
1568 }
1569 
1570 void Hdiv3dPaper201504::ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, std::ostream &out, int DoFT, int DofCond)
1571 {
1572 
1573  int64_t nel = l2mesh->NElements();
1574  int dim = l2mesh->Dimension();
1575  TPZManVector<STATE,10> globalerrors(10,0.);
1576  for (int64_t el=0; el<nel; el++) {
1577  TPZCompEl *cel = l2mesh->ElementVec()[el];
1578  if (!cel) {
1579  continue;
1580  }
1581  TPZGeoEl *gel = cel->Reference();
1582  if (!gel || gel->Dimension() != dim) {
1583  continue;
1584  }
1585  TPZManVector<REAL,10> elerror(10,0.);
1586  elerror.Fill(0.);
1587  cel->EvaluateError(SolExataH1, elerror, 0);
1588 
1589  int nerr = elerror.size();
1590  globalerrors.resize(nerr);
1591 
1592  for (int i=0; i<nerr; i++) {
1593  globalerrors[i] += elerror[i]*elerror[i];
1594  }
1595  }
1596 
1597  out << ndiv << setw(10) << DoFT << setw(20) << DofCond << setw(28) << sqrt(globalerrors[1]) << setw(35) << sqrt(globalerrors[2]) << endl;
1598 
1599 }
1600 
1601 
1602 void Hdiv3dPaper201504::ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
1603 {
1604  int64_t nel = hdivmesh->NElements();
1605  int dim = hdivmesh->Dimension();
1606  TPZManVector<STATE,10> globalerrorsDual(10,0.);
1607  for (int64_t el=0; el<nel; el++) {
1608  TPZCompEl *cel = hdivmesh->ElementVec()[el];
1609  if(cel->Reference()->Dimension()!=dim) continue;
1610  TPZManVector<REAL,10> elerror(10,0.);
1611  elerror.Fill(0.);
1612  cel->EvaluateError(SolExata, elerror, 0);
1613  int nerr = elerror.size();
1614  for (int i=0; i<nerr; i++) {
1615  globalerrorsDual[i] += elerror[i]*elerror[i];
1616  }
1617 
1618 
1619  }
1620 
1621 
1622  nel = l2mesh->NElements();
1623 
1624  TPZManVector<STATE,10> globalerrorsPrimal(10,0.);
1625  for (int64_t el=0; el<nel; el++) {
1626  TPZCompEl *cel = l2mesh->ElementVec()[el];
1627  TPZManVector<REAL,10> elerror(10,0.);
1628  cel->EvaluateError(SolExata, elerror, 0);
1629  int nerr = elerror.size();
1630  globalerrorsPrimal.resize(nerr);
1631 
1632  for (int i=0; i<nerr; i++) {
1633  globalerrorsPrimal[i] += elerror[i]*elerror[i];
1634  }
1635 
1636  }
1637 
1638  // sequence of storage
1639  // for each
1640  errors.PutVal(pos, 0, p); // polinomial order
1641  errors.PutVal(pos, 1, ndiv); // refinement
1642  errors.PutVal(pos, 2, sqrt(globalerrorsPrimal[1])); // pressure
1643  errors.PutVal(pos, 3, sqrt(globalerrorsDual[1])); // flux
1644 
1645 }
1646 
1647 void Hdiv3dPaper201504::ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, std::ostream &out, int DoFT, int DofCond)
1648 {
1649  int64_t nel = hdivmesh->NElements();
1650  int dim = hdivmesh->Dimension();
1651  TPZManVector<STATE,10> globalerrorsDual(10,0.);
1652  for (int64_t el=0; el<nel; el++) {
1653  TPZCompEl *cel = hdivmesh->ElementVec()[el];
1654  if(cel->Reference()->Dimension()!=dim) continue;
1655  TPZManVector<REAL,10> elerror(10,0.);
1656  elerror.Fill(0.);
1657  cel->EvaluateError(SolExata, elerror, 0);
1658  int nerr = elerror.size();
1659  for (int i=0; i<nerr; i++) {
1660  globalerrorsDual[i] += elerror[i]*elerror[i];
1661  }
1662 
1663 
1664  }
1665 
1666 
1667  nel = l2mesh->NElements();
1668 
1669  TPZManVector<STATE,10> globalerrorsPrimal(10,0.);
1670  for (int64_t el=0; el<nel; el++) {
1671  TPZCompEl *cel = l2mesh->ElementVec()[el];
1672  TPZManVector<REAL,10> elerror(10,0.);
1673  cel->EvaluateError(SolExata, elerror, 0);
1674  int nerr = elerror.size();
1675  globalerrorsPrimal.resize(nerr);
1676 
1677  for (int i=0; i<nerr; i++) {
1678  globalerrorsPrimal[i] += elerror[i]*elerror[i];
1679  }
1680 
1681  }
1682 
1683  out << ndiv << setw(10) << DoFT << setw(20) << DofCond << setw(28) << sqrt(globalerrorsPrimal[1]) << setw(35) << sqrt(globalerrorsDual[1]) << endl;
1684 
1685 }
1686 
1687 void Hdiv3dPaper201504::ChangeExternalOrderConnects(TPZCompMesh *mesh){
1688 
1689  int nEl= mesh-> NElements();
1690  int dim = mesh->Dimension();
1691 
1692  int cordermin = -1;
1693  for (int iel=0; iel<nEl; iel++) {
1694  TPZCompEl *cel = mesh->ElementVec()[iel];
1695  if (!cel) continue;
1696  int ncon = cel->NConnects();
1697  int corder = 0;
1698  int nshape = 0;
1699 
1700  if(cel->Dimension()== dim){
1701  for (int icon=0; icon<ncon-1; icon++){
1702  TPZConnect &co = cel->Connect(icon);
1703  corder = co.Order();
1704  nshape = co.NShape();
1705  if(corder!=cordermin){
1706  cordermin = corder-1;
1707  co.SetOrder(cordermin,cel->ConnectIndex(icon));
1708  co.SetNShape(nshape-1);
1709  mesh->Block().Set(co.SequenceNumber(),nshape-1);
1710  }
1711  }
1712  }
1713  }
1714  mesh->ExpandSolution();
1715  mesh->CleanUpUnconnectedNodes();
1716 }
1717 
1718 
1719 
1720 void Hdiv3dPaper201504::SolveSyst(TPZAnalysis &an, TPZCompMesh *fCmesh)
1721 {
1722  std::cout <<"Numero de equacoes "<< fCmesh->NEquations()<< std::endl;
1723 
1724  bool isdirect = true;
1725  bool simetrico = true;
1726  bool isfrontal = false;
1727  if (isdirect)
1728  {
1729  if (simetrico)
1730  {
1731  //TPZSkylineStructMatrix strmat(fCmesh);
1732  if (isfrontal) {
1733  TPZParFrontStructMatrix<TPZFrontSym<STATE> > strmat(fCmesh);
1734  strmat.SetDecomposeType(ELDLt);
1735 
1736  int numthreads = 1;
1737 
1738  strmat.SetNumThreads(numthreads);
1739 
1740  an.SetStructuralMatrix(strmat);
1741  }
1742  else
1743  {
1744  //TPZBandStructMatrix full(fCmesh);
1745  TPZSkylineStructMatrix skylstr(fCmesh); //caso simetrico
1746  // TPZSkylineNSymStructMatrix full(fCmesh);
1747  an.SetStructuralMatrix(skylstr);
1748  }
1749 
1750 
1751  TPZStepSolver<STATE> step;
1752  step.SetDirect(ELDLt); //caso simetrico
1753  an.SetSolver(step);
1754  an.Run();
1755  }
1756  else
1757  {
1758  TPZBandStructMatrix full(fCmesh);
1759  an.SetStructuralMatrix(full);
1760  TPZStepSolver<STATE> step;
1761  step.SetDirect(ELU);
1762  an.SetSolver(step);
1763  an.Run();
1764  }
1765 
1766  }
1767  else
1768  {
1769  TPZSkylineStructMatrix skylstr(fCmesh); //caso simetrico
1770  skylstr.SetNumThreads(10);
1771  an.SetStructuralMatrix(skylstr);
1772 
1773  TPZAutoPointer<TPZMatrix<STATE> > matbeingcopied = skylstr.Create();
1774  TPZAutoPointer<TPZMatrix<STATE> > matClone = matbeingcopied->Clone();
1775 
1776  TPZStepSolver<STATE> *precond = new TPZStepSolver<STATE>(matClone);
1777  TPZStepSolver<STATE> *Solver = new TPZStepSolver<STATE>(matbeingcopied);
1778  precond->SetReferenceMatrix(matbeingcopied);
1779  precond->SetDirect(ELDLt);
1780  Solver->SetGMRES(20, 20, *precond, 1.e-18, 0);
1781  // Solver->SetCG(10, *precond, 1.0e-10, 0);
1782  an.SetSolver(*Solver);
1783  an.Run();
1784  }
1785 
1786 
1787 
1788 }
TPZCompMesh::NElements
int64_t NElements() const
Number of computational elements allocated.
Definition: pzcmesh.h:169
Hdiv3dPaper201504::CMeshMixed
TPZCompMesh * CMeshMixed(TPZGeoMesh *gmesh, TPZVec< TPZCompMesh *> meshvec)
Definition: hdiv3dpaper201504.cpp:1323
TPZCompMesh::Reference
TPZGeoMesh * Reference() const
Returns a pointer to the geometrical mesh associated.
Definition: pzcmesh.h:209
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
ETetraedro
Definition: pzeltype.h:58
TPZCompMesh::SetAllCreateFunctionsContinuous
void SetAllCreateFunctionsContinuous()
Definition: pzcmesh.h:508
TPZConnect::IncrementElConnected
void IncrementElConnected()
Increment fNElConnected.
Definition: pzconnect.h:260
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
TPZBandStructMatrix
Implements Banded Structural Matrices. Structural Matrix.
Definition: pzbstrmatrix.h:18
Hdiv3dPaper201504::MyDoubleComparer
bool MyDoubleComparer(REAL a, REAL b)
Definition: hdiv3dpaper201504.h:199
TPZGeoEl::MaterialId
int MaterialId() const
Returns the material index of the element.
Definition: pzgeoel.h:250
Hdiv3dPaper201504::CreateOneCubo
TPZGeoMesh * CreateOneCubo(int nref)
Definition: hdiv3dpaper201504.cpp:804
TPZGeoMesh::CreateGeoElement
virtual TPZGeoEl * CreateGeoElement(MElementType type, TPZVec< int64_t > &cornerindexes, int matid, int64_t &index, int reftype=1)
Generic method for creating a geometric element. Putting this method centrally facilitates the modifi...
Definition: pzgmesh.cpp:1296
Hdiv3dPaper201504::ErrorH1
void ErrorH1(TPZCompMesh *l2mesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
Definition: hdiv3dpaper201504.cpp:1535
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
Hdiv3dPaper201504::ForcingH1
static void ForcingH1(const TPZVec< REAL > &pt, TPZVec< STATE > &ff, TPZFMatrix< STATE > &flux)
Definition: hdiv3dpaper201504.cpp:1077
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.
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
Hdiv3dPaper201504::tetraedra_2
TPZFMatrix< int > tetraedra_2
Definition: hdiv3dpaper201504.h:130
TPZGeoMesh::NElements
int64_t NElements() const
Number of elements of the mesh.
Definition: pzgmesh.h:129
Hdiv3dPaper201504::CreateOneCuboWithTetraedrons
TPZGeoMesh * CreateOneCuboWithTetraedrons(int64_t nelem)
Definition: hdiv3dpaper201504.cpp:603
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
Hdiv3dPaper201504::CMeshH1
TPZCompMesh * CMeshH1(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdiv3dpaper201504.cpp:1126
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
TPZGeoEl::NodeIndex
virtual int64_t NodeIndex(int i) const =0
Returns the index of the ith node the index is the location of the node in the nodevector of the mesh...
Hdiv3dPaper201504::GMeshWithPrism
TPZGeoMesh * GMeshWithPrism(int ndiv)
Definition: hdiv3dpaper201504.cpp:420
TPZGeoMesh::SetDimension
void SetDimension(int dim)
Set Dimension.
Definition: pzgmesh.h:285
Hdiv3dPaper201504::ErrorPrimalDual
void ErrorPrimalDual(TPZCompMesh *l2mesh, TPZCompMesh *hdivmesh, int p, int ndiv, int pos, TPZFMatrix< REAL > &errors)
Definition: hdiv3dpaper201504.cpp:1602
Hdiv3dPaper201504::Run
void Run(ApproximationSpace problem, Eltype element, TPZVec< int > POrderBeginAndEnd, TPZVec< int > ndivinterval, TPZFMatrix< REAL > &errors)
Definition: hdiv3dpaper201504.cpp:75
TPZCompMesh::LoadReferences
void LoadReferences()
Map this grid in the geometric grid.
Definition: pzcmesh.cpp:482
TPZGeoEl::WhichSide
int WhichSide(TPZVec< int64_t > &SideNodeIds)
Returns the side number which is connected to the SideNodes returns -1 if no side is found...
Definition: pzgeoel.cpp:165
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
Hdiv3dPaper201504::ForcingBC3D
static void ForcingBC3D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1107
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
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
Hdiv3dPaper201504::ForcingBC0D
static void ForcingBC0D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1089
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
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
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
Hdiv3dPaper201504::ChangeExternalOrderConnects
void ChangeExternalOrderConnects(TPZCompMesh *mesh)
Definition: hdiv3dpaper201504.cpp:1687
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
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
TPZCompMesh::SetAllCreateFunctionsDiscontinuous
void SetAllCreateFunctionsDiscontinuous()
Definition: pzcmesh.h:503
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
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
Hdiv3dPaper201504::CMeshPressure
TPZCompMesh * CMeshPressure(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdiv3dpaper201504.cpp:1246
TPZVec::clear
virtual void clear()
Empty the vector, make its size zero.
Definition: pzvec.h:444
TPZCompMesh::SetAllCreateFunctionsHDiv
void SetAllCreateFunctionsHDiv()
Definition: pzcmesh.h:523
substruct_tst14.test.numthreads
list numthreads
Definition: test.py:210
Hdiv3dPaper201504::PrintErrors
void PrintErrors(ApproximationSpace problem, Eltype element, TPZVec< int > POrderBeginAndEnd, TPZVec< int > ndivinterval, TPZVec< REAL > &errors, std::ostream &output)
Definition: hdiv3dpaper201504.cpp:241
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
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
TPZCompMesh::Block
const TPZBlock< STATE > & Block() const
Access the block structure of the solution vector.
Definition: pzcmesh.h:213
Hdiv3dPaper201504::ForcingBC2D
static void ForcingBC2D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1101
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
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.
Hdiv3dPaper201504::SolExataH1
static void SolExataH1(const TPZVec< REAL > &pt, TPZVec< STATE > &solp, TPZFMatrix< STATE > &flux)
Definition: hdiv3dpaper201504.cpp:1056
TPZCompEl::NConnects
virtual int NConnects() const =0
Returns the number of nodes of the element.
TPZGeoElSide::Element
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
TPZGeoElBC
Structure to help the construction of geometric elements along side of a given geometric element...
Definition: pzgeoelbc.h:21
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
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZCompMesh::InsertMaterialObject
int InsertMaterialObject(TPZMaterial *mat)
Insert a material object in the datastructure.
Definition: pzcmesh.cpp:287
Hdiv3dPaper201504::ForcingBC4D
static void ForcingBC4D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1113
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.
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
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
TPZStack
This class implements a stack object. Utility.
Definition: pzcheckmesh.h:14
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
Hdiv3dPaper201504::ForcingBC1D
static void ForcingBC1D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1095
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
Hdiv3dPaper201504::Forcing
static void Forcing(const TPZVec< REAL > &pt, TPZVec< STATE > &ff)
Definition: hdiv3dpaper201504.cpp:1046
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
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
Hdiv3dPaper201504::SolveSyst
void SolveSyst(TPZAnalysis &an, TPZCompMesh *fCmesh)
Definition: hdiv3dpaper201504.cpp:1720
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
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
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
Hdiv3dPaper201504::SolExata
static void SolExata(const TPZVec< REAL > &pt, TPZVec< STATE > &solp, TPZFMatrix< STATE > &flux)
Definition: hdiv3dpaper201504.cpp:1025
Hdiv3dPaper201504::GenerateNodes
void GenerateNodes(TPZGeoMesh *gmesh, int64_t nelem)
Definition: hdiv3dpaper201504.cpp:787
Hdiv3dPaper201504::ForcingBC5D
static void ForcingBC5D(const TPZVec< REAL > &pt, TPZVec< STATE > &disp)
Definition: hdiv3dpaper201504.cpp:1120
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
Hdiv3dPaper201504::CMeshFlux
TPZCompMesh * CMeshFlux(TPZGeoMesh *gmesh, int pOrder, int dim)
Definition: hdiv3dpaper201504.cpp:1187
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