NeoPZ
pzcmesh.cpp
Go to the documentation of this file.
1 
6 #include "pzcmesh.h"
7 #ifdef MACOSX
8 #include <__functional_base> // for less
9 #include <__tree> // for __tree_const_iterator, ope...
10 #endif
11 #include <cmath> // for fabs, sqrt, abs
12 #include <iterator> // for operator!=, reverse_iterator
13 #include <map> // for map, __map_iterator, opera...
14 #include <set> // for set, set<>::reverse_iterator
15 #include <string> // for char_traits, allocator
16 #include <utility> // for pair
17 #include "TPZCompElDisc.h" // for TPZCompElDisc
18 #include "TPZInterfaceEl.h" // for TPZInterfaceElement
19 #ifdef LOG4CXX
20 #include "log4cxx/helpers/objectptr.h" // for ObjectPtrT
21 #include "log4cxx/logger.h" // for Logger
22 #include "log4cxx/propertyconfigurator.h" // for LoggerPtr
23 #endif
24 #include "pzadmchunk.h" // for TPZAdmChunkVector
25 #include "pzblock.h" // for TPZBlock
26 #include "pzbndcond.h" // for TPZBndCond
27 #include "pzcompel.h" // for TPZCompEl, TPZCompElSide
28 #include "pzcondensedcompel.h" // for TPZCondensedCompEl
29 #include "pzconnect.h" // for TPZConnect
30 #include "pzelementgroup.h" // for TPZElementGroup
31 #include "pzeltype.h" // for MElementType::EAgglomerate
32 #include "pzerror.h" // for PZError, DebugStop
33 #include "TPZStream.h" // for TPZStream
34 #include "pzgeoel.h" // for TPZGeoEl
35 #include "pzgeoelside.h" // for TPZGeoElSide
36 #include "pzgmesh.h" // for TPZGeoMesh
37 #include "pzgnode.h" // for TPZGeoNode
38 #include "pzintel.h" // for TPZInterpolatedElement
39 #include "pzinterpolationspace.h" // for TPZInterpolationSpace
40 #include "pzlog.h" // for glogmutex, LOGPZ_DEBUG
41 #include "pzmanvector.h" // for TPZManVector
42 #include "TPZMaterial.h" // for TPZMaterial
43 #include "pzmaterialdata.h" // for TPZSolVec
44 #include "pzmatrix.h" // for TPZFMatrix, TPZMatrix
45 #include "pzmetis.h" // for TPZMetis
46 #include "pzmultiphysicselement.h" // for TPZMultiphysicsElement
47 #include "pzsubcmesh.h" // for TPZSubCompMesh
48 #include "pztransfer.h" // for TPZTransfer
49 #include "pztrnsform.h" // for TPZTransform
50 #include "pzvec.h" // for TPZVec, operator<<
51 
52 #ifndef STATE_COMPLEX
53  #include "TPZAgglomerateEl.h" // for TPZAgglomerateElement
54 #endif
55 
56 #ifdef LOG4CXX
57 static LoggerPtr logger(Logger::getLogger("pz.mesh.tpzcompmesh"));
58 static LoggerPtr aloclogger(Logger::getLogger("pz.allocate"));
59 #endif
60 using namespace std;
61 
62 
64 fElementVec(0),
65 fConnectVec(0),fMaterialVec(),
66 fSolution(0,1) {
67 
68 #ifdef LOG4CXX
69  if (aloclogger->isDebugEnabled()) {
70  std::stringstream sout;
71  sout << "Allocate TPZCompMesh this = " << (void *)this;
72  LOGPZ_DEBUG(aloclogger, sout.str())
73  }
74 #endif
76 
77  //Initializing class members
78  fDimModel = 0;
79  fReference = gr;
80  // fChecked = 0;
81  //fName[0] = '\0';
82  //fName[126] = '\0';
83  if(gr) {
84  SetName( gr->Name() );
85  gr->ResetReference();
86  gr->SetReference(this);
87  SetDimModel(gr->Dimension());
88  }
89  else {
90  SetName( "Computational mesh");
91  }
94 
95  fNmeshes = 0;
96 }
97 
98 
100 fGMesh(gmesh),fElementVec(0),
102 fSolution(0,1)
103 {
104 #ifdef LOG4CXX
105  if (aloclogger->isDebugEnabled()) {
106  std::stringstream sout;
107  sout << "Allocate TPZCompMesh this = " << (void *)this;
108  LOGPZ_DEBUG(aloclogger, sout.str())
109  }
110 #endif
111 
113 
114  //Initializing class members
115  fDimModel = gmesh->Dimension();
116  fReference = gmesh.operator->();
117  // fChecked = 0;
118  //fName[0] = '\0';
119  //fName[126] = '\0';
120  if(fReference) {
121  SetName( fReference->Name() );
123  fReference->SetReference(this);
124  }
125  else {
126  SetName( "Computational mesh");
127  }
130 
131  fNmeshes = 0;
132 }
133 
135 
136 #ifdef LOG4CXX
137  if (aloclogger->isDebugEnabled()) {
138  std::stringstream sout;
139  sout << "Delete TPZCompMesh this = " << (void *)this;
140  LOGPZ_DEBUG(aloclogger, sout.str())
141  }
142 #endif
143 
144 #ifdef LOG4CXX2
145  if (logger->isDebugEnabled()) {
146  std::stringstream sout;
147  Print(sout);
148  LOGPZ_DEBUG(logger, sout.str())
149  }
150 #endif
151  // THIS NEEDS TO INCLUDE THE DELETION ROUTINES OF ALL ITEMS
152  this->CleanUp();
153  TPZGeoMesh * ref = this->Reference();
154  if (ref){
155  if(ref->Reference() == this){
156  ref->ResetReference();
157  }//if == this
158  }//if (ref)
159 }//~
160 
162 
163  // THIS ROUTINE NEEDS TO INCLUDE THE DELETION OF THE LIST POINTERS
164  TPZGeoMesh *ref = Reference();
165  if (ref){
166  ref->ResetReference();
167  this->LoadReferences();
168  }
169 
170  int64_t i, nelem = this->NElements();
171 
172  //deleting subcompmesh
173  for(i=0; i<nelem; i++){
174  TPZCompEl *el = fElementVec[i];
175  TPZSubCompMesh * subc = dynamic_cast<TPZSubCompMesh*>(el);
176  if(subc){
177  delete subc;
178  }
179  }
180 
181 
182  //unwrapping condensed compel
183  for(i=0; i<nelem; i++){
184  TPZCompEl *el = fElementVec[i];
185  TPZCondensedCompEl * cond = dynamic_cast<TPZCondensedCompEl*>(el);
186  if(cond){
187  cond->Unwrap();
188  }
189  }
190 
191  //unwrapping element groups
192  for(i=0; i<nelem; i++){
193  TPZCompEl *el = fElementVec[i];
194  TPZElementGroup * group = dynamic_cast<TPZElementGroup*>(el);
195  if(group){
196  group->Unwrap();
197  }
198  }
199 
200  //deleting interface elements
201  for(i=0; i<nelem; i++){
202  TPZCompEl *el = fElementVec[i];
203  TPZInterfaceElement * face = dynamic_cast<TPZInterfaceElement*>(el);
204  if(face){
205  delete el;
206  }
207  }
208 
209  //deleting other elements
210  for(i=0; i<nelem; i++) {
211  TPZCompEl *el = fElementVec[i];
212  if(!el) continue;
213  delete el;
214  }
215 
216  fElementVec.Resize(0);
218  fConnectVec.Resize(0);
220  nelem = NMaterials();
221  std::map<int, TPZMaterial *>::iterator it;
222  for (it = fMaterialVec.begin(); it != fMaterialVec.end(); it++) {
223  delete it->second;
224  }
225  fMaterialVec.clear();
226 
227  fBlock.SetNBlocks(0);
229  fSolution.Redim(0,1);
230 }
231 
232 void TPZCompMesh::SetName (const string &nm) {
233  fName = nm;
234 }
235 
236 void TPZCompMesh::Print (std::ostream & out) const {
237 
238  //ComputeNodElCon();
239  out << "\n\t\tCOMPUTABLE GRID INFORMATIONS:\n\n";
240  out << "TITLE-> " << fName << "\n\n";
241 
242  out << "number of connects = " << NConnects() << std::endl;
243  out << "number of elements = " << NElements() << std::endl;
244  out << "number of materials = " << NMaterials() << std::endl;
245  out << "dimension of the mesh = " << this->Dimension() << std::endl;
246 
247  out << "\n\t Connect Information:\n\n";
248  int64_t i, nelem = NConnects();
249  for(i=0; i<nelem; i++) {
250  if(fConnectVec[i].SequenceNumber() == -1) {
251  if(fConnectVec[i].HasDependency()) {
252  cout << " TPZCompMesh::Print inconsistency of connect\n";
253  cout << " Index " << i << ' ';
254  fConnectVec[i].Print(*this,std::cout);
255  }
256  continue;
257  }
258  out << " Index " << i << ' ';
259  fConnectVec[i].Print(*this,out);
260  }
261  out << "\n\t Computable Element Information:\n\n";
262  nelem = NElements();
263  for(i=0; i<nelem; i++) {
264  if(!fElementVec[i]) continue;
265  TPZCompEl *el = fElementVec[i];
266  out << "\n Index " << i << ' ';
267  el->Print(out);
268  TPZMultiphysicsElement *mpel = dynamic_cast<TPZMultiphysicsElement *>(el);
269  if(!mpel){
270  if(!el->Reference()) continue;
271  out << "\tReference Index = " << el->Reference()->Index() << std::endl << std::endl;
272  }
273  }
274  out << "\n\t Material Information:\n\n";
275  std::map<int, TPZMaterial * >::const_iterator mit;
276  nelem = NMaterials();
277  for(mit=fMaterialVec.begin(); mit!= fMaterialVec.end(); mit++) {
278  TPZMaterial *mat = mit->second;
279  if (!mat) {
280  DebugStop();
281  }
282  mat->Print(out);
283  }
284 }
285 
288  if(!mat) return -1;
289  int matid = mat->Id();
290  if (fMaterialVec.find(matid) != fMaterialVec.end()) {
291  DebugStop();
292  }
293  fMaterialVec[matid] = mat;
294  return fMaterialVec.size();
295 }
296 
297 TPZMaterial * TPZCompMesh::FindMaterial(int matid){ // find the material object with id matid
298  TPZMaterial * result = 0;
299  std::map<int, TPZMaterial * >::iterator mit;
300  mit = fMaterialVec.find(matid);
301  if(mit != fMaterialVec.end())
302  {
303  result = mit->second;
304  }
305  return result;
306 }
307 
308 void TPZCompMesh::AutoBuild(const std::set<int> *MaterialIDs) {
309 #ifdef PZDEBUG
310  {
311  TPZGeoMesh *gmesh = Reference();
312  TPZCheckGeom check(gmesh);
313  check.CheckIds();
314  }
315 #endif
316  if(MaterialIDs)
317  {
318  fCreate.BuildMesh(*this, *MaterialIDs);
319  }
320  else {
321  fCreate.BuildMesh(*this);
322  }
323 
324 }
325 
326 void TPZCompMesh::AutoBuildContDisc(const TPZVec<TPZGeoEl*> &continuous, const TPZVec<TPZGeoEl*> &discontinuous) {
327 #ifdef PZDEBUG
328  {
329  TPZGeoMesh *gmesh = Reference();
330  TPZCheckGeom check(gmesh);
331  check.CheckIds();
332  }
333 #endif
334 
336  int64_t nelem = elvec.NElements();
337 
338  int64_t neltocreate = 0;
339  int64_t index;
340  for(int64_t i=0; i<nelem; i++) {
341  TPZGeoEl *gel = elvec[i];
342  if(!gel) continue;
343  if(!gel->HasSubElement()) {
344  neltocreate++;
345  }
346  }
347 
348  int64_t nbl = fBlock.NBlocks();
349  if(neltocreate > nbl) fBlock.SetNBlocks(neltocreate);
350  fBlock.SetNBlocks(nbl);
351 
352  //Creating continuous elements
354  int64_t ncont = continuous.NElements();
355  for(int64_t i = 0; i < ncont; i++){
356  TPZGeoEl *gel = continuous[i];
357  if(!gel) continue;
358  if(!gel->HasSubElement()) {
359  int printing = 0;
360  if (printing) {
361  gel->Print(cout);
362  }
363 
364  if(gel->NumInterfaces() == 0){
365  CreateCompEl(gel,index);
366  }
367  }
368  }
369 
370  //Creating discontinuous elements
372  int64_t ndisc = discontinuous.NElements();
373  for(int64_t i = 0; i < ndisc; i++){
374  TPZGeoEl *gel = discontinuous[i];
375  if(!gel) continue;
376  if(!gel->HasSubElement()) {
377  int printing = 0;
378  if (printing) {
379  gel->Print(cout);
380  }
381 
382  if(gel->NumInterfaces() == 0){
383  CreateCompEl(gel,index);
384  }
385  }
386  }
387 
388  this->InitializeBlock();
389 }
390 
392  ExpandSolution();
394 }
395 
397  fBlock.Resequence();
398  int64_t ibl,nblocks = fBlock.NBlocks();
399 
400  //TPZFMatrix<REAL> OldSolution(fSolution);
401  TPZFMatrix<STATE> OldSolution(fSolution);
402 
403  int64_t cols = fSolution.Cols();
404  fSolution.Redim(fBlock.Dim(),cols);
405  int64_t minblocks = nblocks < fSolutionBlock.NBlocks() ? nblocks : fSolutionBlock.NBlocks();
406  /*
407  int ic;
408  for(ic=0; ic<cols; ic++) {
409  for(ibl = 0;ibl<minblocks;ibl++) {
410  int oldsize = fSolutionBlock.Size(ibl);
411  int oldposition = fSolutionBlock.Position(ibl);
412  int newsize = fBlock.Size(ibl);
413  int newposition = fBlock.Position(ibl);
414  int minsize = (oldsize < newsize) ? oldsize : newsize;
415  int ieq;
416  int offset = 0;
417  if(Discontinuous)offset = newsize - oldsize;
418  for(ieq=0; ieq<minsize; ieq++) {
419  fSolution(newposition+ieq+offset,ic) = OldSolution(oldposition+ieq,ic);
420  }
421  }
422  }
423  */
424  int64_t ic;
425  for(ic=0; ic<cols; ic++) {
426  for(ibl = 0;ibl<minblocks;ibl++) {
427  int64_t oldsize = fSolutionBlock.Size(ibl);
428  int64_t oldposition = fSolutionBlock.Position(ibl);
429  int64_t newsize = fBlock.Size(ibl);
430  int64_t newposition = fBlock.Position(ibl);
431  int64_t minsize = (oldsize < newsize) ? oldsize : newsize;
432  int64_t ieq;
433  for(ieq=0; ieq<minsize; ieq++) {
434  fSolution(newposition+ieq,ic) = OldSolution(oldposition+ieq,ic);
435  }
436  }
437  }
439 }
440 
442 
443  int64_t nrow = mat.Rows();
444  int64_t ncol = mat.Cols();
445  int64_t solrow = fSolution.Rows();
446  fSolution.Resize(solrow, ncol);
447  int64_t i,j;
448  STATE val;
449  for(j=0;j<ncol;j++)
450  {
451  for(i=0;i<nrow;i++)
452  {
453  val = (mat.GetVal(i,j));
454  fSolution(i,j) = val;
455  }
456 
457  }
458 
459 
460  int64_t nelem = NElements();
461  TPZCompEl *cel;
462  for(i=0; i<nelem; i++) {
463  cel = fElementVec[i];
464  if(!cel) continue;
465  cel->LoadSolution();
466  }
467 
468 }
469 
471 {
472  int64_t nel = this->NElements();
473  for (int64_t iel = 0; iel < nel; iel++) {
474  TPZCompEl *cel = this->Element(iel);
475  if (!cel) {
476  continue;
477  }
479  }
480 }
481 
483 
484  // Reference()->ResetReference();
485  Reference()->SetReference(this);
486  int64_t i, nelem = NElements();
487  for(i=0; i<nelem; i++) {
488  TPZCompEl *el = fElementVec[i];
489  if(!el) continue;
490  el->LoadElementReference();
491  /* TPZGeoEl *gel = el->Reference();
492  if(!gel) continue;
493  gel->SetReference(el);
494  */
495  }
496 }
497 
499  ComputeNodElCon();
500  int64_t i, nconnects = NConnects();
501  int64_t ndepblocks = 0, nvalidblocks = 0, nremoved = 0, ncondensed = 0;
502  for (i=0;i<nconnects;i++)
503  {
504  TPZConnect &no = fConnectVec[i];
505  int64_t seq = no.SequenceNumber();
506  if(!no.NElConnected() && seq != -1)
507  {
508  nremoved++;
509  }
510  else if(!no.HasDependency() && no.NElConnected() && !no.IsCondensed()) nvalidblocks++;
511  else if(!no.HasDependency() && no.NElConnected() && no.IsCondensed()) ncondensed++;
512  else if(no.HasDependency() && no.NElConnected()) ndepblocks++;
513  }
514  int need = 0;
515  for (i=0;i<nconnects;i++) {
516  TPZConnect &no = fConnectVec[i];
517  if (no.SequenceNumber() == -1) continue;
518  if (no.HasDependency() && no.NElConnected() == 0) {
519  PZError << "TPZCompMesh::CleanUpUnconnectedNodes node has dependency\n";
520  continue;
521  }
522  if (!no.NElConnected() && no.SequenceNumber() != -1)
523  {
524  need = 1;
525  break;
526  }
527  if (no.HasDependency() && no.SequenceNumber() < nvalidblocks)
528  {
529  need = 1;
530  break;
531  } else if(!no.HasDependency() && !no.IsCondensed() && no.SequenceNumber() >= nvalidblocks)
532  {
533  need = 1;
534  break;
535  } else if(!no.HasDependency() && no.IsCondensed() && no.SequenceNumber() < nvalidblocks)
536  {
537  need = 1;
538  break;
539  }
540  }
541  int64_t nblocks = fBlock.NBlocks();
542  TPZManVector<int64_t> permute(nblocks,-1), down(nblocks,0);
543  int64_t idepblocks = 0, iremovedblocks= 0, icondensed = 0;
544 
545  if (need) {
546  for(i=0; i<nconnects; i++) {
547  TPZConnect &no = fConnectVec[i];
548  if(no.SequenceNumber() == -1) continue;
549  int seq = no.SequenceNumber();
550  if(no.NElConnected() == 0)
551  {
552  permute[seq] = nvalidblocks+ndepblocks+iremovedblocks+ncondensed;
553  down[seq] = 1;
554  fBlock.Set(seq,0);
555  no.Reset();
556  // no.SetSequenceNumber(-1);
557  fConnectVec.SetFree(i);
558  iremovedblocks++;
559  }
560  else if(no.HasDependency()) {
561  permute[seq] = nvalidblocks+ncondensed+idepblocks;
562  down[seq] = 1;
563  idepblocks++;
564  }
565  else if(no.IsCondensed())
566  {
567  permute[seq] = nvalidblocks+icondensed;
568  down[seq] = 1;
569  icondensed++;
570 
571  }
572  }
573  for(i=1; i<nblocks; i++) down[i] += down[i-1];
574  for(i=0; i<nblocks; i++)
575  {
576  if(permute[i] == -1)
577  {
578  permute[i] = i-down[i];
579  }
580  }
581  }
582 #ifdef LOG4CXX
583  if(need)
584  if (logger->isDebugEnabled())
585  {
586  std::stringstream sout;
587  sout << "permute to put the free connects to the back\n";
588  if(nblocks < 50)
589  {
590  sout << "original sequence numbers|nelconected\n";
591  int64_t nel = fConnectVec.NElements();
592  for (int64_t el=0; el<nel; el++) {
593  TPZConnect &c = fConnectVec[el];
594  int64_t seqnum = c.SequenceNumber();
595  sout << seqnum << '|' << c.NElConnected() << " ";
596  }
597  sout << std::endl;
598  }
599  if(nblocks < 50) {
600  for (i=0;i<nblocks;i++) sout << permute[i] << ' ';
601  sout << std::endl;
602  }
603  sout << "need = " << need << endl;
604  LOGPZ_DEBUG(logger,sout.str());
605  }
606 #endif
607 
608  if (need) {
609 #ifdef PZDEBUG
610  std::set<int64_t> check;
611  nconnects = permute.NElements();
612  for(i=0; i<nconnects; i++) check.insert(permute[i]);
613  if(static_cast<int>(check.size()) != nconnects)
614  {
615  cout << __PRETTY_FUNCTION__ << " The permutation vector is not a permutation!\n" << permute << endl;
616  DebugStop();
617  }
618 #endif
619  Permute(permute);
620 
621 #ifdef LOG4CXX
622  if (logger->isDebugEnabled() && nblocks < 50)
623  {
624  if(nblocks < 50)
625  {
626  std::stringstream sout;
627  sout << "after permute sequence numbers|nelconected\n";
628  int64_t nel = fConnectVec.NElements();
629  for (int64_t el=0; el<nel; el++) {
630  TPZConnect &c = fConnectVec[el];
631  int64_t seqnum = c.SequenceNumber();
632  sout << seqnum << '|' << c.NElConnected() << " ";
633  }
634  LOGPZ_DEBUG(logger, sout.str())
635  }
636 
637  }
638 #endif
639  int64_t nel = fConnectVec.NElements();
640  for (int64_t i=0;i<nel;i++) {
641  TPZConnect &no = fConnectVec[i];
642  if (no.NElConnected() == 0 && no.SequenceNumber() >= nblocks-nremoved) {
643  no.Reset();
644  fConnectVec.SetFree(i);
645  }
646  else if(no.NElConnected() == 0 && no.SequenceNumber() != -1)
647  {
648  DebugStop();
649  }
650  }
651 #ifdef PZDEBUG
652  {
653  int64_t nel = fConnectVec.NElements();
654  for (int64_t el=0; el<nel; el++) {
655  TPZConnect &c = fConnectVec[el];
656  int64_t seqnum = c.SequenceNumber();
657  if (seqnum > nblocks-nremoved) {
658  DebugStop();
659  }
660  }
661  }
662 #endif
663 
664  fBlock.SetNBlocks(nblocks-nremoved);
665  }
666 }
667 
669 
670  int64_t i, nelem = NConnects();
671  for(i=0; i<nelem; i++) {
672  TPZConnect &no = fConnectVec[i];
673  if(no.SequenceNumber() == -1) continue;
674  no.ResetElConnected();
675  }
676 
677 
678  TPZStack<int64_t> nodelist;
679  int64_t numnod;
680  // modified Philippe 22/7/97
681  // in order to account for constrained nodes
682  nelem = NElements();
683  for(i=0; i<nelem; i++) {
684  TPZCompEl *el = fElementVec[i];
685  if(!el) continue;
686  nodelist.Resize(0);
687  el->BuildConnectList(nodelist);
688  numnod = nodelist.NElements();
689  for (int64_t in=0; in<numnod; ++in) {
690  int64_t dfnindex = nodelist[in];
691  TPZConnect *dfn = &fConnectVec[dfnindex];
692  dfn->IncrementElConnected();
693  }
694  }
695 }
696 
697 void TPZCompMesh::ComputeNodElCon(TPZVec<int> &nelconnected ) const {
698 
699  int64_t i, nelem = NConnects();
700  nelconnected.Resize(nelem);
701  nelconnected.Fill(0);
702  TPZStack<int64_t> nodelist;
703  int64_t numnod;
704  // modified Philippe 22/7/97
705  // in order to account for constrained nodes
706  nelem = NElements();
707  for(i=0; i<nelem; i++) {
708  TPZCompEl *el = fElementVec[i];
709  if(!el) continue;
710  nodelist.Resize(0);
711  el->BuildConnectList(nodelist);
712  numnod = nodelist.NElements();
713  for (int64_t in=0; in<numnod; ++in) {
714  int64_t dfnindex = nodelist[in];
715  nelconnected[dfnindex]++;
716  }
717  }
718 }
719 
720 
722  int64_t neq = 0;
723  int64_t i, ncon = NConnects();
724  for(i=0; i<ncon; i++) {
725  TPZConnect &df = fConnectVec[i];
726  if(df.HasDependency() || df.IsCondensed() || !df.NElConnected() || df.SequenceNumber() == -1){
727  continue;
728  }
729 
730  int dofsize = df.NShape()*df.NState();
731 #ifdef PZDEBUG
732  // check the consistency between the block size and the data structure of the connect
733  {
734  int64_t seqnum = df.SequenceNumber();
735  int64_t blsize = fBlock.Size(seqnum);
736  if (blsize != dofsize) {
737  DebugStop();
738  }
739  }
740 #endif
741  neq += dofsize;
742  }
743  return neq;
744 }
745 
748 
749  int bw = 0;
750  TPZStack<int64_t> connectlist;
751  // modified Philippe 24/7/97
752  // in order to take dependent nodes into account
753 
754  int64_t i, nelem = NElements();
755  for(i=0; i<nelem; i++) {
756  connectlist.Resize(0);
757  TPZCompEl *el = fElementVec[i];
758  if(!el) continue;
759  el->BuildConnectList(connectlist);
760  int64_t nnod = connectlist.NElements();
761  if(!nnod) continue;
762  // look for a node which has equations associated with it
763  int64_t ifirstnode = 0;
764  TPZConnect *np = &fConnectVec[connectlist[ifirstnode++]];
765  while(ifirstnode < nnod && (np->HasDependency() || np->IsCondensed() || !fBlock.Size(np->SequenceNumber()))) {
766  np = &fConnectVec[connectlist[ifirstnode++]];
767  }
768  int64_t ibl = np->SequenceNumber();
769  int64_t loweq = fBlock.Position(ibl);
770  int64_t higheq = loweq+fBlock.Size(ibl)-1;
771  for(int64_t n=ifirstnode;n<nnod;n++) {
772  np = &fConnectVec[connectlist[n]];
773  if(np->HasDependency() || np->IsCondensed() ) continue;
774  int64_t ibl = np->SequenceNumber();
775  if(!fBlock.Size(ibl)) continue;
776  int64_t leq = fBlock.Position(ibl);
777  int64_t heq = leq+fBlock.Size(ibl)-1;
778  loweq = (loweq > leq) ? leq : loweq;
779  higheq = (higheq < heq) ? heq : higheq;
780  }
781  int64_t elbw = higheq - loweq;
782  bw = (bw < elbw) ? elbw : bw;
783  }
784  return bw;
785 }
786 
788 
789  TPZStack<int64_t> connectlist;
790  // modified Philippe 24/7/97
791  // in order to take dependent nodes into account
792 
793 #ifdef PZDEBUG
794  TPZAdmChunkVector <TPZConnect > & connectVec = this->ConnectVec();
795  int maxSequenceNumberIndependentConnect = 0;
796  TPZVec<int> depConInd(0,0);//list of connects that have dependencies
797  for (int i = 0; i < connectVec.NElements(); i++) {
798  if (connectVec[i].HasDependency() ){
799  int oldSize = depConInd.size();
800  depConInd.Resize( oldSize + 1 );
801  depConInd[oldSize] = i;
802  continue;
803  }
804  if (connectVec[i].SequenceNumber() > maxSequenceNumberIndependentConnect && !connectVec[i].IsCondensed())
805  {
806  maxSequenceNumberIndependentConnect = connectVec[i].SequenceNumber();
807  }
808  }
809  for (int i = 0; i < depConInd.size(); i++) {
810  if (connectVec[ depConInd[i] ].SequenceNumber() < maxSequenceNumberIndependentConnect ) {
811 
812  std::cout<<"A connect that has dependency has "
813  <<"a sequency number smaller than an independent connect. "
814  <<"have you tried this->CleanUpUnconnectedNodes() after "
815  <<"building the computational mesh?"<<std::endl;
816  DebugStop();
817  }
818  }
819 #endif
820 
821  int64_t neq = NEquations();
822  skyline.Resize(neq);
823  if (neq == 0) {
824  return;
825  }
826 // cout << "Element skyline\n";
827  //int eleq=0;
828  int64_t i, n, l, nelem = NElements();
829  for(i=0; i<neq; i++) skyline[i] = i;
830  for(i=0; i<nelem; i++) {
831  TPZCompEl *el = fElementVec[i];
832  if(!el) continue;
833  // if(!el) continue;
834  connectlist.Resize(0);
835  el->BuildConnectList(connectlist);
836  int64_t nnod = connectlist.NElements();
837  if(!nnod) continue;
838  // look for a connect with global equations associated to it
839  int64_t ifirstnode = 0;
840  TPZConnect *np = &fConnectVec[connectlist[0]];
841  while(ifirstnode < nnod && (np->HasDependency() || np->IsCondensed()) ) {
842  ifirstnode++;
843  if (ifirstnode == nnod) {
844  break;
845  }
846  np = &fConnectVec[connectlist[ifirstnode]];
847  }
848  int64_t ibl = np->SequenceNumber();
849  int64_t loweq = fBlock.Position(ibl);
850  int64_t higheq = loweq+fBlock.Size(ibl)-1;
851  for(n=ifirstnode;n<nnod;n++) {
852  np = &fConnectVec[connectlist[n]];
853  if(np->HasDependency() || np->IsCondensed()) continue;
854  int64_t ibl = np->SequenceNumber();
855  int64_t leq = fBlock.Position(ibl);
856  int64_t heq = leq+fBlock.Size(ibl)-1;
857  //for(int _eq=leq; _eq<= heq; _eq++) {
858  // if((eleq%20==0)) cout << endl;
859  // cout << _eq << ' ';
860  // eleq++;
861  //}
862  loweq = (loweq > leq) ? leq : loweq;
863  higheq = (higheq < heq) ? heq : higheq;
864  }
865 // std::cout << "Element " << i << " loweq " << loweq << " higheq " << higheq << std::endl;
866 // cout << "Equations ";
867  for(n=ifirstnode;n<nnod;n++) {
868  np = &fConnectVec[connectlist[n]];
869  if(np->HasDependency() || np->IsCondensed()) continue;
870  int64_t ibl = np->SequenceNumber();
871  int64_t leq = fBlock.Position(ibl);
872  int64_t heq = leq+fBlock.Size(ibl);
873  for(l=leq;l<heq;l++) {
874  skyline[l] = skyline[l] < loweq ? skyline[l] : loweq;
875 // cout << l << "/" << skyline[l] << " ";
876  }
877  }
878 // cout << endl;
879  }
880 // std::cout << "Skyline " << skyline << std::endl;
881 }
882 
884 
885  //TPZBlock<REAL> &localblock = Block();
886  TPZBlock<STATE> &localblock = Block();
887  //TPZBlock<REAL> &coarseblock = coarsemesh.Block();
888  TPZBlock<STATE> &coarseblock = coarsemesh.Block();
889  // adapt the block size of the blocks, dividing by the number of variables
890  // of the material
891  int nmat = NMaterials();
892  if(!nmat) {
893  PZError << "TPZCompMesh::BuildTransferMatrix, no material object found\n";
894  return;
895  }
896  TPZMaterial * mat;
897  mat = fMaterialVec.begin()->second;
898  int nvar = mat->NStateVariables();
899  int dim = mat->Dimension();
900 
901  transfer.SetBlocks(localblock,coarseblock,nvar,NIndependentConnects(),coarsemesh.NIndependentConnects());
903  coarsemesh.LoadReferences();
904  int64_t nelem = NElements();
905  int64_t i;
906  for(i=0; i<nelem; i++) {
907  if(!fElementVec[i]) continue;
908  TPZInterpolationSpace * finecel = dynamic_cast<TPZInterpolationSpace *> (fElementVec[i]);
909  if(!finecel) continue;
910  if(finecel->Dimension() != dim) continue;
911  TPZGeoEl *finegel = finecel->Reference();
912  TPZGeoEl *coarsegel = finegel;
913  if(!finegel) {
914  cout << "TPZCompMesh::BuildTransferMatrix is not implemented for super elements\n";
915  continue;
916  }
917 
918  while(coarsegel && !coarsegel->Reference()) {
919  coarsegel = coarsegel->Father();
920  }
921  if(!coarsegel) {
922  cout << "TPZCompMesh::BuildTransferMatrix corresponding coarse element not found\n";
923  finecel->Print(cout);
924  continue;
925  }
926 
927  TPZInterpolationSpace * coarsel = dynamic_cast<TPZInterpolationSpace *> ( coarsegel->Reference() );
928  if(!coarsel) continue;
929 
930  if(coarsel->Mesh() != &coarsemesh) {
931  cout << "TPZCompMesh::BuildTransferMatrix is not implemented for transfers"
932  " between superelements\n";
933  continue;
934  }
935  TPZTransform<> t(coarsel->Dimension());
936  t=finegel->BuildTransform2(finegel->NSides()-1,coarsegel,t);
937  finecel->BuildTransferMatrix(*coarsel,t,transfer);
938  }
939 }
940 
942  TPZTransfer<STATE> &transfer) {
943 #ifdef STATE_COMPLEX
944  DebugStop();
945 #else
946  //TPZBlock<REAL> &localblock = Block();
947  TPZBlock<STATE> &localblock = Block();
948  //TPZBlock<REAL> &transferblock = transfermesh.Block();
949  TPZBlock<STATE> &transferblock = transfermesh.Block();
950  // adapt the block size of the blocks, dividing by the number of variables
951  // of the material
952  int nmat = NMaterials();
953  if(!nmat) {
954  PZError << "TPZCompMesh::BuildTransferMatrixDesc no material object found\n";
955  return;
956  }
957  TPZMaterial * mat;
958  mat = fMaterialVec.begin()->second;
959  int nvar = mat->NStateVariables();
960  int dim = mat->Dimension();
961  //o seguinte �igual ao nmero de conects da malha
962  int64_t ncon = NIndependentConnects(),coarncon = transfermesh.NIndependentConnects();
963  transfer.SetBlocks(localblock,transferblock,nvar,ncon,coarncon);
964  Reference()->ResetReference();//geom�ricos apontam para nulo
965  transfermesh.LoadReferences();
966  //geom�ricos apontam para computacionais da malha atual
967  TPZAgglomerateElement *aggel = 0;
968  TPZAdmChunkVector<TPZCompEl *> &elvec = transfermesh.ElementVec();
969  int64_t nelem = elvec.NElements();
970  int64_t i;
971  for(i=0; i<nelem; i++) {
972  TPZCompEl *comp = elvec[i];
973  if(!comp) continue;
974  if(comp->Dimension() != dim) continue;
975  if(comp->Type() != EAgglomerate){
976  PZError << "TPZCompMesh::BuildTransferMatrixDesc mesh agglomerated"
977  << " with element of volume not agglomerated\n";
978  continue;
979  }
980  aggel = dynamic_cast<TPZAgglomerateElement *>(comp);
981  //TPZStack<int> elvec;
982  //retorna todos os descont�uos aglomerados por aggel
983  //aggel->IndexesDiscSubEls(elvec);
984  //int size = elvec.NElements(),i;
985  int64_t size = aggel->NIndexes(),i;
986  for(i=0;i<size;i++){
987  //TPZCompElDisc *disc = dynamic_cast<TPZCompElDisc *>(fElementVec[elvec[i]]);
988  // TPZCompElDisc *disc = dynamic_cast<TPZCompElDisc *>(aggel->FineElement(i));
989  TPZCompElDisc *disc = dynamic_cast<TPZCompElDisc *>(aggel->SubElement(i));
990  if(!disc){
991  PZError << "TPZCompMesh::BuildTransferMatrixDesc index with null"
992  << " elemento\n";
993  continue;
994  }
995  if(disc->Type() != EDiscontinuous) {
996  PZError << "TPZCompMesh::BuildTransferMatrixDesc index of not"
997  << " discontinous element\n";
998  continue;
999  }
1000  disc->BuildTransferMatrix(*aggel,transfer);
1001  }
1002  }
1003 #endif
1004 }
1005 
1006 /*
1007  void TPZCompMesh::CreateConnectBC() {
1008  TPZGeoMesh *geo = Reference();
1009  TPZAdmChunkVector<TPZGeoNodeBC> *geobndcondvec = &geo->BCNodeVec();
1010  int ibc, nnodebc = geobndcondvec->NElements();
1011  for(ibc=0; ibc<nnodebc; ibc++) {
1012  TPZGeoNodeBC &gbc = (*geobndcondvec)[ibc];
1013  int bcnumber = gbc.fBCId;
1014  TPZMaterial *bc = FindMaterial(bcnumber);
1015  if(!bc) continue;
1016  // find a geometric element which has a
1017  // corresponding computational element
1018  TPZGeoElSide gel(gbc.fGeoEl,gbc.fGeoElSide);
1019  TPZStack<TPZGeoElSide> neighbourset;
1020  gel.AllNeighbours(neighbourset);
1021  int in=0,nneigh;
1022  nneigh = neighbourset.NElements();
1023  while(in<nneigh && !neighbourset[in].Reference().Exists()) in++;
1024  // neighbour = gel.Neighbour();
1025  // while(neighbour.Exists() && neighbour != gel && !neighbour.Reference().Exists()) {
1026  // neighbour = neighbour.Neighbour();
1027  // }
1028  // if(!neighbour.Exists() || ! neighbour.Reference().Exists()) continue;
1029  if(in == nneigh) continue;
1030  TPZCompElSide cel = neighbourset[in].Reference();
1031  TPZConnect *df = &cel.Element()->Connect(neighbourset[in].Side());
1032  if(!df) continue;
1033  TPZConnectBC dfbc(df,(TPZBndCond *) bc);
1034  int key = BCConnectVec().AllocateNewElement();
1035  BCConnectVec()[key] = dfbc;
1036  }
1037  }
1038  */
1039 
1040 /*
1041  void TPZCompMesh::ComputeConnecttoElGraph(TPZVec<int> &firstel, TPZVec<int> &connectelgraph){
1042  int connectstackstore[50];
1043  TPZStack<int> connectstack(connectstackstore,50);
1044  int i, ncon = NConnects();
1045  firstel.Resize(ncon+1);
1046  firstel[0] = 0;
1047  for(i=0; i<ncon; i++) {
1048  TPZConnect &c = fConnectVec[i];
1049  int seqnum = c.SequenceNumber();
1050  if(seqnum == -1) {
1051  firstel[i+1] = firstel[i];
1052  } else {
1053  firstel[i+1] = firstel[i] + c.NElConnected();
1054  }
1055  }
1056  connectelgraph.Resize(firstel[ncon]);
1057  connectelgraph.Fill(-1);
1058  int nelem = NElements();
1059  for(i=0; i<nelem; i++) {
1060  TPZCompEl *el = fElementVec[i];
1061  if(!el) continue;
1062  connectstack.Resize(0);
1063  el->BuildConnectList(connectstack);
1064  int in;
1065  ncon = connectstack.NElements();
1066  for(in=0; in<ncon; in++) {
1067  int ic = connectstack[in];
1068  int first = firstel[ic];
1069  int last = firstel[ic+1];
1070  while(connectelgraph[first] != -1 && first < last) first++;
1071  if(first == last) {
1072  PZError << "TPZCompMesh::ComputeConnecttoElGraph wrong data structure\n";
1073  continue;
1074  }
1075  connectelgraph[first] = i;
1076  }
1077  }
1078  }
1079  */
1081  int64_t i, ncon = NConnects();
1082  int64_t NIndependentConnects = 0;
1083  for(i=0; i<ncon; i++) {
1084  TPZConnect &c = fConnectVec[i];
1085  if(c.HasDependency() || c.IsCondensed() || c.SequenceNumber() == -1) continue;
1086  NIndependentConnects++;
1087  }
1088  return NIndependentConnects;
1089 }
1090 
1092 
1093  std::set<int> mat_ids;
1094  for (auto item : fMaterialVec) {
1095  mat_ids.insert(item.first);
1096  }
1097  ComputeElGraph(elgraph, elgraphindex, mat_ids);
1098 }
1099 
1100 void TPZCompMesh::ComputeElGraph(TPZStack<int64_t> &elgraph, TPZVec<int64_t> &elgraphindex, std::set<int> & mat_ids){
1101 
1102  int64_t i, ncon;
1103  TPZStack<int64_t> connectstack;
1104  int64_t nelem = NElements();
1105  elgraphindex.Resize(nelem+1);
1106  elgraphindex[0] = 0;
1107  elgraph.Resize(0);
1108  int64_t curel=0;
1109  int64_t nindep = this->NIndependentConnects();
1110  for(i=0; i<nelem; i++) {
1111  TPZCompEl *el = fElementVec[i];
1113  bool has_material_Q = true;
1114  if (el) {
1115  has_material_Q = el->HasMaterial(mat_ids);
1116  }
1117  if(!el || !has_material_Q){
1118  elgraphindex[curel+1]=elgraph.NElements();
1119  curel++;
1120  continue;
1121  }
1122 
1123 
1124 
1125  connectstack.Resize(0);
1126  el->BuildConnectList(connectstack);
1127  int64_t in;
1128  ncon = connectstack.NElements();
1129  for(in=0; in<ncon; in++) {
1130  int ic = connectstack[in];
1131  TPZConnect &c = fConnectVec[ic];
1132  if(c.HasDependency() || c.IsCondensed()) continue;
1133 #ifdef PZDEBUG
1134  if (c.SequenceNumber() >= nindep) {
1135  DebugStop();
1136  }
1137 #endif
1138  elgraph.Push(c.SequenceNumber());
1139  }
1140  elgraphindex[curel+1]=elgraph.NElements();
1141  curel++;
1142  }
1143  elgraphindex.Resize(curel+1);
1144 }
1145 
1146 void TPZCompMesh::Divide(int64_t index,TPZVec<int64_t> &subindex,int interpolate) {
1147 
1148  TPZCompEl * el = fElementVec[index];
1149  if (!el) {
1150  PZError << "TPZCompMesh::Divide element not found index = " << index << endl;
1151  subindex.Resize(0);
1152  return;
1153  }
1154  if(index != el->Index()){
1155  PZError << "TPZCompMesh::Divide - element->Index() != index " << endl;
1156  subindex.Resize(0);
1157  return;
1158  }
1159 
1160 
1161  el->Divide(index,subindex,interpolate);
1162 }
1163 
1164 void TPZCompMesh::Coarsen(TPZVec<int64_t> &elements, int64_t &index, bool CreateDiscontinuous) {
1165  int64_t i;
1166  const int64_t nelem = elements.NElements();
1167 
1168  if(!nelem) {
1169  index = -1;
1170  return;
1171  }//if
1172 
1173  TPZGeoEl *father = 0;
1174  TPZCompEl *cel = fElementVec[elements[0]];
1175  if (!cel) {
1176  index = -1;
1177  return;
1178  }//if
1179 
1180  if(cel) father = cel->Reference()->Father();
1181  if(!father) {
1182  index = -1;
1183  return;
1184  }//if
1185 
1186  for(i=1;i<nelem;i++) {
1187  if(!fElementVec[elements[i]]) {
1188  index = -1;
1189  return;
1190  }//if
1191 
1192  TPZGeoEl *father2 = fElementVec[elements[i]]->Reference()->Father();
1193  if(!father2 || father != father2) {
1194  index = -1;
1195  return;
1196  }//if
1197 
1198  }//for i
1199 
1200  if(nelem != father->NSubElements()){
1201  cout << "TPZCompEl::Coarsen : incomplete list of elements sons\n";
1202  }//if
1203 
1204  for(i=0; i<nelem; i++) {
1205  TPZInterpolationSpace * cel = dynamic_cast<TPZInterpolationSpace*>(this->ElementVec()[elements[i]]);
1206  if (!cel) continue;
1207  cel->RemoveInterfaces();
1208  TPZInterpolatedElement * intel = dynamic_cast<TPZInterpolatedElement *>(cel);
1209  if (intel){
1211  }//if (intel)
1212 
1213  cel->Reference()->ResetReference();
1214  }//for
1215 
1216 
1217  for(i=0; i<nelem; i++) {
1218  TPZCompEl * cel = this->ElementVec()[elements[i]];
1219  if (cel) delete cel;
1220  }
1221 
1222  if (CreateDiscontinuous) fCreate.SetAllCreateFunctionsDiscontinuous();
1224 
1225  TPZCompEl * newcel = CreateCompEl(father,index);
1226 
1227  TPZCompElDisc * newdisc = dynamic_cast<TPZCompElDisc*>(newcel);
1228  if (newdisc){
1229  newdisc->SetDegree( this->GetDefaultOrder() );
1230  }
1231 
1232 }//method
1233 
1234 void TPZCompMesh::Discontinuous2Continuous(int64_t disc_index, int64_t &new_index) {
1235 
1236  TPZInterpolationSpace *cel = dynamic_cast<TPZInterpolationSpace*> (fElementVec[disc_index]);
1237  if (!cel) {
1238  new_index = -1;
1239  return;
1240  }//if
1241 
1242  TPZCompElDisc * disc = dynamic_cast<TPZCompElDisc*>(cel);
1243  if (!disc) {
1244  new_index = -1;
1245  return;
1246  }//if
1247 
1248  TPZGeoEl * ref = cel->Reference();
1249  if (!ref){
1250  LOGPZ_FATAL(logger, "Computational element without reference");
1251  return;
1252  }
1253  cel->RemoveInterfaces();
1254  cel->Reference()->ResetReference();
1255  // this->fElementVec[ cel->Index() ] = NULL;
1256  // delete cel;
1257 
1259  TPZCompEl * newcel = CreateCompEl(ref,new_index);
1260  TPZInterpolatedElement * intel = dynamic_cast<TPZInterpolatedElement*>(newcel);
1261  intel->CreateInterfaces(false);
1262 
1263  if (!intel){
1264  LOGPZ_FATAL(logger, "New created element is not an interpolated element as requested.");
1265  }
1266 
1267  if(!ref->Reference()){
1268  LOGPZ_FATAL(logger, "New created element is not referenced by geometric element");
1269  }
1270 
1271  if (ref->Reference() != newcel){
1272  LOGPZ_FATAL(logger, "New created element is not the same as referenced by geometric element");
1273  }
1274 
1275  ExpandSolution();
1276  intel->InterpolateSolution(*disc);
1277  delete cel;
1278 
1279 }//method
1280 
1282 
1283  int64_t n = this->ElementVec().NElements();
1284 
1285  for(int64_t i = 0; i < n; i++){
1286  TPZCompElDisc * disc = dynamic_cast<TPZCompElDisc*>( this->ElementVec()[i] );
1287  if (!disc) continue;
1288  disc->RemoveInterfaces();
1289  }//for
1290 
1291 #ifdef PZDEBUG
1292  {
1293  n = this->ElementVec().NElements();
1294  for(int64_t i = 0; i < n; i++){
1295  TPZCompEl * cel = this->ElementVec()[i];
1296  if (!cel) continue;
1297  // MElementType type = cel->Type();
1298  TPZInterfaceElement * face = dynamic_cast<TPZInterfaceElement*>(cel);
1299  if(face){
1300  PZError << __PRETTY_FUNCTION__ << " - At this point no TPZInterfaceElement may exist.\n";
1301  }
1302  }
1303  }
1304 #endif
1305 
1306  n = this->ElementVec().NElements();
1307  for(int64_t i = 0; i < n; i++){
1308  TPZCompElDisc * disc = dynamic_cast<TPZCompElDisc*>( this->ElementVec()[i] );
1309  if (!disc) continue;
1310  disc->CreateInterfaces();
1311  }//for
1312 
1313 }//method
1314 
1315 #ifdef LOG4CXX
1316 #include "pzsubcmesh.h"
1317 #endif
1318 
1320 // it is a gather permutation
1322 
1323  ExpandSolution();
1324  // if (permute.NElements() != fBlock.NBlocks()) {
1325  // PZError << "TPZCompMesh::Permute : permute vector size not equal to fBlock size\n";
1326  // }
1327 #ifdef LOG4CXX
1328  if(logger->isDebugEnabled())
1329  {
1330  std::stringstream sout;
1331  TPZSubCompMesh *submesh = dynamic_cast<TPZSubCompMesh *> (this);
1332  if (submesh) {
1333  sout << "Index = " << submesh->Index() << " ";
1334  }
1335  sout << "Permutation " << permute;
1336  std::set<int64_t> permset;
1337  if (permute.size() != 0) {
1338  permset.insert(&permute[0],(&permute[permute.size()-1]+1));
1339  }
1340  sout << " permute size " << permute.size() << " no dist elements " << permset.size();
1341  LOGPZ_DEBUG(logger,sout.str())
1342  }
1343 #endif
1344  int64_t i,j;
1345  int64_t permutenel = permute.NElements();
1346  for (i = 0; i < permutenel; i++) fBlock.Set(permute[i],fSolutionBlock.Size(i));
1347  fBlock.Resequence();
1348  if (fSolution.Rows() != 0) {
1349  //TPZFMatrix<REAL> newsol(fSolution);
1350  TPZFMatrix<STATE> newsol(fSolution);
1351  for (i=0;i<fBlock.NBlocks();i++) {
1352  int64_t oldpos = fSolutionBlock.Position(i);
1353  int64_t newpos;
1354  if(i < permutenel) {
1355  newpos = fBlock.Position(permute[i]);
1356  } else {
1357  newpos = fBlock.Position(i);
1358  }
1359  for (j=0;j<fSolutionBlock.Size(i);j++) fSolution(newpos+j,0) = newsol(oldpos+j,0);
1360  } //a sol. inicial esta em newsol
1361  }
1362 
1364  int64_t ncon = NConnects();
1365  for(i=0; i<ncon; i++) {
1366  TPZConnect &df = fConnectVec[i];
1367  int64_t seqnum = df.SequenceNumber();
1368  if(seqnum == -1) continue;
1369  if(seqnum < permutenel) df.SetSequenceNumber(permute[seqnum]);
1370  }
1371 }
1372 
1373 void TPZCompMesh::ConnectSolution(std::ostream & out) {
1374 
1375  out << "\n\t\tCONNECT INDEX SOLUTION:\n\n";
1376  int64_t i;
1377  int64_t ncon = NConnects();
1378  for(i=0; i<ncon; i++) {
1379  out << i << ") ";
1380  TPZConnect &df = ConnectVec()[i];
1381  int64_t seqnum = df.SequenceNumber();
1382  if(df.NElConnected()==0) {
1383  out << "free node" << endl;
1384  } else if (seqnum < 0 || Block().Size(seqnum)==0) {
1385  out << "non solution connect" << endl;
1386  } else {
1387  int64_t pos = Block().Position(seqnum);
1388  for(int64_t j=0;j<Block().Size(seqnum);j++)
1389  out << Solution()(pos+j,0) << " ";
1390  out << std::endl;
1391  }
1392  }
1393 }
1394 
1395 void TPZCompMesh::EvaluateError(std::function<void (const TPZVec<REAL> &loc,TPZVec<STATE> &val,TPZFMatrix<STATE> &deriv)> fp, bool store_error, TPZVec<REAL> &errorSum) {
1396 
1397  errorSum.Resize(3);
1398  errorSum.Fill(0.);
1399 
1400  TPZManVector<REAL,3> true_error(3);
1401  true_error.Fill(0.);
1402 
1403 
1404  TPZCompEl *cel;
1405 
1406  //soma de erros sobre os elementos
1407  for(int64_t el=0;el< fElementVec.NElements();el++) {
1408  cel = fElementVec[el];
1409  if(!cel || cel->Material()->Id() < 0) continue;
1410  cel->EvaluateError(fp,true_error,store_error);
1411 
1412  int64_t nerrors = true_error.NElements();
1413  errorSum.Resize(nerrors,0.);
1414  for(int64_t ii = 0; ii < nerrors; ii++)
1415  errorSum[ii] += true_error[ii]*true_error[ii];
1416  }
1417 
1418  int64_t nerrors = errorSum.NElements();
1419  for(int64_t ii = 0; ii < nerrors; ii++)
1420  errorSum[ii] = sqrt(errorSum[ii]);
1421 }
1422 
1424  int changed = 1;
1425  while(changed) {
1426  changed = 0;
1427  int64_t nel = fElementVec.NElements();
1428  int64_t el;
1429  TPZVec<int64_t> subelindex;
1430  for(el=0; el<nel; el++) {
1432  TPZCompEl *elp = fElementVec[el];
1433 
1434  if(!elp || !dynamic_cast<TPZInterpolatedElement*>(elp) ) continue;
1435 
1436  TPZMaterial * mat = elp->Material();
1437  // this statement determines thata the element is associated with a boundary condition
1438  TPZBndCond *bnd = dynamic_cast<TPZBndCond *>(mat);
1439  if(!bnd) continue;
1440  // if(mat && mat->Id() >= 0) continue;
1441  int nsides = elp->Reference()->NSides();
1442  int is;
1443  int nc = elp->Reference()->NCornerNodes();
1444  for(is=nc; is<nsides; is++) {
1445  TPZCompElSide elpside(elp,is);
1446  // this should never be called
1447  if(elp->Reference()->SideDimension(is) == 0) continue;
1448  elvec.Resize(0);
1449  // verify if there are smaller elements connected to this boundary element
1450  elpside.HigherLevelElementList(elvec,0,0);
1451  if(elvec.NElements()) {
1452  // the first small element
1453  TPZGeoElSide fatherside = elvec[0].Reference();//el BC
1454  // elp is the element we will eventually refine
1455  int face = elp->Reference()->NSides() - 1;
1456  // this is the volume side of the element
1457  while(fatherside.Exists()) {
1458  if(elp->Reference()->NeighbourExists(face,fatherside)) break;
1459  fatherside = fatherside.Father2();
1460  }
1461  // fatherside is a neighbour of the current element
1462  // I wouldnt know when this test could fail
1463  if(fatherside.Exists()) {
1464 #ifdef LOG4CXX
1465  if(logger->isDebugEnabled())
1466  {
1467  std::stringstream sout;
1468  sout << "Dividing element " << el << " of type " << elp->Reference()->TypeName();
1469  LOGPZ_DEBUG(logger,sout.str().c_str());
1470  }
1471 #endif
1472  Divide(el,subelindex,0);
1473  changed = 1;
1474  break;
1475  }
1476  }
1477  // we are working on the last side and nothing was divided
1478  if(is == nsides-1) {
1479  TPZInterpolatedElement *elpint = dynamic_cast <TPZInterpolatedElement *> (elp);
1480  if(!elpint) continue;
1481  int porder = elpint->PreferredSideOrder(is);
1482  int maxorder = 0;
1483  elvec.Resize(0);
1484  elpside.EqualLevelElementList(elvec,0,0);
1485  int64_t eq;
1486  for(eq=0; eq<elvec.NElements(); eq++) {
1487  TPZInterpolatedElement *eqel = dynamic_cast<TPZInterpolatedElement *> (elvec[eq].Element());
1488  int eqside = elvec[eq].Side();
1489  if(!eqel) continue;
1490  if(maxorder < eqel->PreferredSideOrder(eqside)) maxorder = eqel->PreferredSideOrder(eqside);
1491  }
1492  // set the order to the largest order of all connecting elements
1493  if(porder < maxorder) {
1494 #ifdef LOG4CXX
1495  if(logger->isDebugEnabled())
1496  {
1497  std::stringstream sout;
1498  sout << "Refining element " << el << " to order " << maxorder;
1499  LOGPZ_DEBUG(logger,sout.str().c_str());
1500  }
1501 #endif
1502  elpint->PRefine(maxorder);
1503  changed = 1;
1504  }
1505  }
1506  }
1507  }
1508  }
1509  InitializeBlock();
1510 }
1511 
1512 int64_t TPZCompMesh::PutinSuperMesh (int64_t local, TPZCompMesh *super){
1513  if (super != this) return -1;
1514  else return local;
1515 }
1516 
1517 int64_t TPZCompMesh::GetFromSuperMesh (int64_t superind, TPZCompMesh *super){
1518  if (super != this) return -1;
1519  else return superind;
1520 }
1521 
1522 REAL TPZCompMesh::CompareMesh(int var, char *matname){
1523 
1524  REAL error = 0.;
1525  int64_t i=0;
1526  for (i=0;i<fElementVec.NElements();i++){
1527  TPZCompEl *el = fElementVec[i];
1528  if(el) error+= el->CompareElement(var,matname);
1529  }
1530  return (error);
1531 }
1532 
1534  if(sol.NElements() != NElements()) {
1535  cout << "TPZCompMesh::SetElementSolution size of the vector doesn't match\n";
1536  }
1537 
1538 #ifdef LOG4CXX
1539  if(logger->isDebugEnabled())
1540  {
1541  std::stringstream sout;
1542  STATE norm=0.;
1543  for (int64_t ii=0; ii<sol.size(); ii++) {
1544  norm += sol[ii];
1545  }
1546  norm = sqrt(norm);
1547  sout << "Norma da solucao " << i << " norma " << norm;
1548  LOGPZ_DEBUG(logger, sout.str())
1549  }
1550 #endif
1552  int64_t el,nel= NElements();
1553  for(el=0; el<nel; el++) {
1554  fElementSolution(el,i) = sol[el];
1555  }
1556 }
1557 
1558 void TPZCompMesh::GetRefPatches(std::set<TPZGeoEl *> &grpatch){
1559  int64_t i;
1560  int64_t nel = NElements();
1561  // cout << "GetRefPatches\n" << nel << endl;
1562  for (i=0; i<nel; i++){
1563  if (fElementVec[i]){
1564  TPZGeoEl *gel = fElementVec[i]->Reference();
1565  if (gel) gel = fElementVec[i]->GetRefElPatch();
1566  if (gel)
1567  {
1568  grpatch.insert(gel);
1569  }
1570  }
1571  }
1572  return;
1573 }
1574 
1575 
1576 void TPZCompMesh::GetNodeToElGraph(TPZVec<int64_t> &nodtoelgraph, TPZVec<int64_t> &nodtoelgraphindex, TPZStack<int64_t> &elgraph, TPZVec<int64_t> &elgraphindex){
1577 
1578  ComputeElGraph(elgraph,elgraphindex);
1579 
1580  // int i,j;
1581  // for (i=0; i<elgraphindex.NElements()-1; i++){
1582  // cout << "Block: " << i << "\t";
1583  // for (j = elgraphindex[i]; j<elgraphindex[i+1]; j++){
1584  // cout << elgraph[j] << "\t";
1585  // }
1586  // cout << endl;
1587  // }
1588 
1589  TPZMetis renum(elgraphindex.NElements() -1 ,NIndependentConnects());
1590  renum.SetElementGraph(elgraph,elgraphindex);
1591  renum.NodeToElGraph(elgraph,elgraphindex,nodtoelgraph, nodtoelgraphindex);
1592  /* TPZRenumbering *re = &renum; */
1593  /* re->Print(nodtoelgraph, nodtoelgraphindex , "Grapho de n� para elementos "); */
1594 
1595 
1596 }
1597 
1598 
1599 void TPZCompMesh::GetElementPatch(TPZVec<int64_t> nodtoelgraph, TPZVec<int64_t> nodtoelgraphindex, TPZStack<int64_t> &elgraph, TPZVec<int64_t> &elgraphindex,int64_t elind ,TPZStack<int64_t> &patch){
1600 
1601  // int aux =0;
1602  //TPZAVLMap<int,int> elconmap(aux);
1603  std::set<int64_t > elconmap;
1604  int64_t i,j;
1605  for (i= elgraphindex[elind]; i<elgraphindex[elind+1];i++){
1606  int64_t node = elgraph[i];
1607  for (j = nodtoelgraphindex[node]; j<nodtoelgraphindex[node+1]; j++){
1608  elconmap.insert(nodtoelgraph[j]);
1609 
1610  }
1611  }
1612  patch.Resize(0);
1613 
1614  //TPZPix iter = elconmap.First();
1615  set<int64_t >::iterator iter = elconmap.begin();
1616  while(iter!=elconmap.end()){
1617  //patch.Push(elconmap.Key(iter));
1618  patch.Push(*iter);//elconmap.Key(iter));
1619  iter++;//elconmap.Next(iter);
1620  }
1621 }
1622 
1627 fSolution(copy.fSolution), fBlock(copy.fBlock),
1629 {
1630 #ifdef LOG4CXX
1631  if (aloclogger->isDebugEnabled()) {
1632  std::stringstream sout;
1633  sout << "Allocate TPZCompMesh this = " << (void *)this;
1634  LOGPZ_DEBUG(aloclogger, sout.str())
1635  }
1636 #endif
1637 
1642  copy.CopyMaterials(*this);
1643  int64_t nel = copy.fElementVec.NElements();
1644  fElementVec.Resize(nel);
1645  int64_t iel;
1646  for(iel = 0; iel<nel; iel++) fElementVec[iel] = 0;
1647  for(iel = 0; iel<nel; iel++) {
1648  TPZCompEl *cel = copy.fElementVec[iel];
1649  if(cel && !dynamic_cast<TPZInterfaceElement* >(cel) )
1650  {
1651  /*TPZCompEl *clone = */ cel->Clone(*this);
1652  /*#ifdef LOG4CXX
1653  {
1654  std::stringstream sout;
1655  sout << "original\n";
1656  cel->Print(sout);
1657  sout << "cloned\n";
1658  clone->Print(sout);
1659  LOGPZ_DEBUG(logger,sout.str())
1660  }
1661  #endif*/
1662  }
1663  }
1665  ComputeNodElCon();
1666 // int nconn = fConnectVec.NElements();
1667 // for(iel=0;
1668  /*#ifdef LOG4CXX
1669  {
1670  std::stringstream sout;
1671  Print(sout);
1672  LOGPZ_DEBUG(logger,sout.str())
1673  }
1674  #endif*/
1675 // for(iel = 0; iel<nel; iel++) {
1676 // TPZCompEl *cel = copy.fElementVec[iel];
1677 // if(cel && dynamic_cast<TPZInterfaceElement* >(cel) ) cel->Clone(*this);
1678 // }
1679  fDimModel = copy.fDimModel;
1680  fName = copy.fName;
1681 }
1682 
1684 {
1685  CleanUp();
1686  fReference = copy.fReference;
1688  fConnectVec = copy.fConnectVec;
1689  copy.CopyMaterials(*this);
1691  fSolution = copy.fSolution;
1693  fBlock = copy.fBlock;
1697  int64_t nel = copy.fElementVec.NElements();
1698  fElementVec.Resize(nel);
1699  int64_t iel;
1700  for(iel = 0; iel<nel; iel++) fElementVec[iel] = nullptr;
1701  for(iel = 0; iel<nel; iel++) {
1702  TPZCompEl *cel = copy.fElementVec[iel];
1703  if(cel && !dynamic_cast<TPZInterfaceElement* >(cel) )
1704  {
1705  cel->Clone(*this);
1706  }
1707  }
1708  for(iel = 0; iel<nel; iel++)
1709  {
1710  TPZCompEl *cel = copy.fElementVec[iel];
1711  if(cel && dynamic_cast<TPZInterfaceElement* >(cel) ) cel->Clone(*this);
1712  }
1713  fDimModel = copy.fDimModel;
1714  fName = copy.fName;
1715  fCreate = copy.fCreate;
1716 
1717  return *this;
1718 }
1719 
1721  return new TPZCompMesh(*this);
1722 }
1723 /*
1724  TPZGeoMesh *gmesh = Reference();
1725  gmesh->ResetReference();
1726  TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
1727  CopyMaterials(cmesh);
1728 
1729  int iel, nel = fElementVec.NElements();
1730  for (iel=0;iel<nel;iel++){
1731  if (!fElementVec[iel]) continue;
1732  TPZInterpolatedElement *intel = dynamic_cast<TPZInterpolatedElement *> (fElementVec[iel]);
1733 
1734  TPZGeoEl *gel = fElementVec[iel]->Reference();
1735  if (!gel) continue;
1736  int indexcomp;
1737  gel->CreateCompEl(*cmesh,indexcomp);
1738 
1739  TPZInterpolatedElement *clintel = dynamic_cast<TPZInterpolatedElement *> (cmesh->ElementVec()[indexcomp]);
1740  int side, nsides = intel->NConnects();
1741  for (side =0;side<nsides;side++){
1742  int porder = intel->SideOrder(side);
1743  clintel->PRefine(side,porder);
1744  int clorder = clintel->SideOrder(side);
1745  if(porder != clorder) {
1746  cout << "PZCompMesh::Clone order was not honoured side " <<
1747  side << " porder " << porder << " clorder " << clorder <<
1748  endl;
1749  }
1750  int elndof = intel->Connect(side).NDof(*this);
1751  int clelndof = clintel->Connect(side).NDof(*cmesh);
1752  if(elndof != clelndof) {
1753  cout << "PZCompMesh::Clone ndof different side " <<
1754  side << " porder " << porder << " clorder " << clorder <<
1755  " elndof = "<< elndof << " clelndof " << clelndof << endl;
1756  }
1757  }
1758  }
1759 
1760  cmesh->InitializeBlock();
1761 
1762  gmesh->ResetReference();
1763  LoadReferences();
1764 
1765 
1766  int clnel = cmesh->ElementVec().NElements();
1767  // if (clnel != nel){
1768  // cout << "Clone mesh inconsistancy: clone elements failed!\n";
1769  // return 0;
1770  // }
1771 
1772  for (iel=0;iel<clnel;iel++){
1773  TPZCompEl *cel = cmesh->ElementVec()[iel];
1774  TPZGeoEl *gel = cel->Reference();
1775  TPZCompEl *el = gel->Reference();
1776  int ic,nc = el->NConnects();
1777  for (ic=0;ic<nc;ic++){
1778  int elseqnum = el->Connect(ic).SequenceNumber();
1779  int elndof = el->Connect(ic).NDof(*this);
1780  int clelseqnum = cel->Connect(ic).SequenceNumber();
1781  int clelndof = cel->Connect(ic).NDof(*cmesh);
1782  if (elndof != clelndof){
1783  cout << "Number of degree of freedom incompatible between clone and original mesh!\n";
1784  cout << "Clone connect id: " << ic << " Clone dof: " << clelndof << " Original dof: " << clelndof << endl;
1785  continue;
1786  }
1787  int idof;
1788  for (idof=0; idof<elndof; idof++){
1789  cmesh->fSolutionBlock(clelseqnum,0,idof,0) = Block()(elseqnum,0,idof,0);
1790  }
1791  }
1792  }
1793  return cmesh;
1794  }
1795  */
1796 
1798  // int nmat = fMaterialVec.size();
1799  std::map<int, TPZMaterial * >::const_iterator mit;
1800  // int m;
1801  for(mit=fMaterialVec.begin(); mit!=fMaterialVec.end(); mit++) {
1802  if(!dynamic_cast<TPZBndCond *> (mit->second)) mit->second->Clone(mesh.fMaterialVec);
1803  }
1804  for(mit=fMaterialVec.begin(); mit!=fMaterialVec.end(); mit++) {
1805  if(dynamic_cast<TPZBndCond *> (mit->second)) mit->second->Clone(mesh.fMaterialVec);
1806  }
1807 
1808 }
1809 
1811 
1812  int64_t nel = ElementVec().NElements(),i,j;
1813  if(nel == 0) cout << "\nTPZCompMesh::DeltaX no elements\n";
1814  REAL maxdist = 0.0,dist=0.0;
1815  TPZVec<REAL> point0(3),point1(3);
1816  TPZGeoNode *node0,*node1;
1817  for(i=0;i<nel;i++){
1818  TPZCompEl *com = ElementVec()[i];
1819  if(!com) continue;
1820  if(com->Type() == EInterface || com->Material()->Id() < 0) continue;
1821  node0 = com->Reference()->NodePtr(0);
1822  node1 = com->Reference()->NodePtr(1);
1823  for(j=0;j<3;j++){
1824  point0[j] = node0->Coord(j);
1825  point1[j] = node1->Coord(j);
1826  }
1827  dist = TPZGeoEl::Distance(point0,point1);
1828  if(dist > maxdist) maxdist = dist;
1829  }
1830  return maxdist;
1831 }
1832 
1834 
1835  int64_t nel = ElementVec().NElements(),i;
1836  if(nel == 0) cout << "\nTPZCompMesh::MaximumRadiusOfMesh no elements\n";
1837  REAL maxdist = 0.0,dist=0.0;
1838  TPZVec<REAL> point0(3),point1(3);
1839  for(i=0;i<nel;i++){
1840  TPZCompEl *com = ElementVec()[i];
1841  if(!com) continue;
1842  if(com->Type() == EInterface || com->Material()->Id() < 0) continue;
1843  dist = com->MaximumRadiusOfEl();
1844  if(dist > maxdist) maxdist = dist;
1845  }
1846  return maxdist;
1847 }
1848 
1850 
1851  int64_t nel = ElementVec().NElements(),i;
1852  if(nel == 0) cout << "\nTPZCompMesh::MaximumRadiusOfMesh no elements\n";
1853  REAL mindist =10000.0,dist=0.0;
1854  for(i=0;i<nel;i++){
1855  TPZCompEl *com = ElementVec()[i];
1856  if(!com) continue;
1857  int type = com->Type();
1858  if( type == EInterface || com->Material()->Id() < 0 ) continue;
1859  if(type == EAgglomerate) dist = dynamic_cast<TPZCompElDisc *>(com)->LesserEdgeOfEl();
1860  else dist = com->LesserEdgeOfEl();
1861  if(dist < mindist) mindist = dist;
1862  }
1863  return mindist;
1864 }
1865 
1869 void TPZCompMesh::ComputeFillIn(int64_t resolution, TPZFMatrix<REAL> &fillin){
1870  ComputeNodElCon();
1871  int64_t nequations = NEquations();
1872  int64_t divider = nequations/resolution;
1873  if(divider*resolution != nequations) divider++;
1874  REAL factor = 1./(divider*divider);
1875  fillin.Redim(resolution,resolution);
1876 
1877  TPZStack<int64_t> graphelindex, graphel, graphnodeindex, graphnode;
1878  this->ComputeElGraph(graphel,graphelindex);
1880  renum.ConvertGraph(graphel,graphelindex,graphnode,graphnodeindex);
1881  std::map<int64_t,TPZConnect *> seqtoconnect;
1882  int ic,ncon = fConnectVec.NElements();
1883  for(ic=0; ic<ncon; ic++) {
1884  TPZConnect &c = fConnectVec[ic];
1885  if(c.HasDependency() || c.IsCondensed() || c.SequenceNumber() < 0) continue;
1886  seqtoconnect[c.SequenceNumber()] = &c;
1887  }
1888  int64_t iseqnum;
1889  for(iseqnum = 0; iseqnum < graphnodeindex.NElements()-1; iseqnum++) {
1890  if(!seqtoconnect.count(iseqnum)) continue;
1891  int64_t firstieq = Block().Position(iseqnum);
1892  int64_t lastieq = Block().Size(iseqnum)+firstieq;
1893  int64_t firstnode = graphnodeindex[iseqnum];
1894  int64_t lastnode = graphnodeindex[iseqnum+1];
1895  {
1896  int64_t ieq;
1897  for(ieq=firstieq; ieq<lastieq; ieq++) {
1898  int64_t rowp = ieq/divider;
1899  int64_t ieq2;
1900  for(ieq2=firstieq; ieq2<lastieq; ieq2++) {
1901  int64_t rowp2 = ieq2/divider;
1902  fillin(rowp,rowp2) += factor;
1903  }
1904  }
1905  }
1906  int64_t in;
1907  for(in=firstnode; in<lastnode; in++) {
1908  int64_t jseqnum = graphnode[in];
1909  int64_t firstjeq = Block().Position(jseqnum);
1910  int64_t lastjeq = Block().Size(jseqnum)+firstjeq;
1911  int64_t ieq;
1912  for(ieq=firstieq; ieq<lastieq; ieq++) {
1913  int64_t rowp = ieq/divider;
1914  int64_t jeq;
1915  for(jeq=firstjeq; jeq<lastjeq; jeq++) {
1916  int64_t colp = jeq/divider;
1917  fillin(rowp,colp) += factor;
1918  }
1919  }
1920  }
1921  }
1922 }
1924 
1925  // * * A MALHA ATUAL DEVE SER AGLOMERADA * * *
1926 
1927  // TPZBlock &localblock = Block();
1928  // TPZBlock &transferblock = finemesh.Block();
1929  // adapt the block size of the blocks, dividing by the number of variables
1930  // of the material
1931  int64_t neq = NEquations();
1932  projectsol.Redim(neq,1);
1933  projectsol.Zero();
1934  int nmat = NMaterials();
1935  if(!nmat) {
1936  PZError << "TPZCompMesh::BuildTransferMatrixDesc2 no material object found\n";
1937  return;
1938  }
1939  Reference()->ResetReference();//geom�ricos apontam para nulo
1940  LoadReferences();
1941 
1942 #ifdef STATE_COMPLEX
1943  DebugStop();
1944 #else
1945  TPZMaterial * mat = fMaterialVec.begin()->second;
1946  //geom�ricos apontam para computacionais da malha atual
1947  int dim = mat->Dimension();
1948  TPZAgglomerateElement *aggel = 0;
1950  int64_t nelem = elvec.NElements();
1951 
1952  for(int64_t i=0; i<nelem; i++) {
1953  TPZCompEl *comp = elvec[i];
1954  if(!comp) continue;
1955  if(comp->Dimension() != dim) continue;
1956  if(comp->Type() != EAgglomerate){
1957  PZError << "TPZCompMesh::BuildTransferMatrixDesc2 mesh agglomerated"
1958  << " with element of volume not agglomerated\n";
1959  continue;
1960  }
1961  aggel = dynamic_cast<TPZAgglomerateElement *>(comp);
1962  aggel->ProjectSolution(projectsol);
1963  }
1964 #endif
1965 }
1970  return Hash("TPZCompMesh");
1971 }
1975 void TPZCompMesh::Write(TPZStream &buf, int withclassid) const { //ok
1977  TPZPersistenceManager::WritePointer(fGMesh.operator->(), &buf);
1978  buf.Write(&fName);
1980  fConnectVec.Write(buf, withclassid);
1981  std::map<int, TPZMaterial*> internal_materials;
1982  std::map<int, TPZMaterial*> boundary_materials;
1983  for (auto mat_pair : fMaterialVec) {
1984  if (dynamic_cast<TPZBndCond*>(mat_pair.second)){
1985  boundary_materials.insert(mat_pair);
1986  } else {
1987  internal_materials.insert(mat_pair);
1988  }
1989  }
1990  buf.WritePointers(internal_materials);
1991  buf.WritePointers(boundary_materials);
1992  fSolutionBlock.Write(buf,0);
1993  fSolution.Write(buf,0);
1994  fBlock.Write(buf,0);
1995  fElementSolution.Write(buf,0);
1996  buf.Write(&fDimModel);
1997  buf.Write(&fDefaultOrder);
1998  fCreate.Write(buf, withclassid);
1999  buf.Write(&fNmeshes);
2000 
2001 }
2002 
2006 void TPZCompMesh::Read(TPZStream &buf, void *context) { //ok
2007  fReference = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::GetInstance(&buf));
2008  fGMesh = TPZAutoPointerDynamicCast<TPZGeoMesh >(TPZPersistenceManager::GetAutoPointer(&buf));
2009  buf.Read(&fName);
2011  fConnectVec.Read(buf, context);
2012  buf.ReadPointers(fMaterialVec); //internal materials
2013  buf.ReadPointers(fMaterialVec); //boundary materials
2014  fSolutionBlock.Read(buf, NULL);
2015  fSolution.Read(buf,NULL);
2016  fBlock.Read(buf, NULL);
2017  fElementSolution.Read(buf, NULL);
2018  buf.Read(&fDimModel);
2019  buf.Read(&fDefaultOrder);
2020  fCreate.Read(buf, context);
2021  buf.Read(&fNmeshes);
2022 }
2023 
2024 #include "TPZGeoElement.h"
2025 #include "pzgeoelrefless.h"
2026 #include "TPZGeoCube.h"
2027 #include "pzshapecube.h"
2028 #include "TPZRefCube.h"
2029 #include "pzshapelinear.h"
2030 #include "TPZGeoLinear.h"
2031 #include "TPZRefLinear.h"
2032 #include "pzrefquad.h"
2033 #include "pzshapequad.h"
2034 #include "pzgeoquad.h"
2035 #include "pzshapetriang.h"
2036 #include "pzreftriangle.h"
2037 #include "pzgeotriangle.h"
2038 #include "pzshapeprism.h"
2039 #include "pzrefprism.h"
2040 #include "pzgeoprism.h"
2041 #include "pzshapetetra.h"
2042 #include "pzreftetrahedra.h"
2043 #include "pzgeotetrahedra.h"
2044 #include "pzshapepiram.h"
2045 #include "pzrefpyram.h"
2046 #include "pzgeopyramid.h"
2047 #include "pzrefpoint.h"
2048 #include "pzgeopoint.h"
2049 #include "pzshapepoint.h"
2050 
2051 
2052 void TPZCompMesh::ConvertDiscontinuous2Continuous(REAL eps, int opt, int dim, TPZVec<STATE> &celJumps){
2053  const int64_t nelements = this->NElements();
2054  celJumps.Resize(nelements);
2055  int64_t numbersol = Solution().Cols();
2056  TPZVec<TPZCompEl*> AllCels(nelements);
2057  AllCels.Fill(NULL);
2058  celJumps.Fill(0.0);
2059 
2060  for(int64_t i = 0; i < nelements; i++){
2061  TPZCompEl * cel = this->ElementVec()[i];
2062  if (!cel) continue;
2063  TPZInterfaceElement * face = dynamic_cast<TPZInterfaceElement *>(cel);
2064  if (!face){
2065  AllCels[i] = cel;
2066  continue;
2067  }
2068  TPZSolVec facejump;
2069  face->EvaluateInterfaceJump(facejump,opt);
2070  const int64_t leftel = face->LeftElement()->Index();
2071  const int64_t rightel = face->RightElement()->Index();
2072 #ifdef PZDEBUG
2073  if(this->ElementVec()[leftel] != face->LeftElement()){
2074  LOGPZ_FATAL(logger, "inconsistent data structure");
2075  DebugStop();
2076  }
2077  if(this->ElementVec()[rightel] != face->RightElement()){
2078  LOGPZ_FATAL(logger, "inconsistent data structure");
2079  DebugStop();
2080  }
2081 #endif
2082 
2083  STATE jumpNorm = 0.;
2084  for (int64_t is=0; is<numbersol; is++) {
2085  for(int64_t ij = 0; ij < facejump.NElements(); ij++){
2086  jumpNorm += facejump[is][ij]*facejump[is][ij];
2087  }
2088  }
2089  jumpNorm = sqrt(jumpNorm);
2090 
2091  celJumps[leftel] += jumpNorm;
2092  celJumps[rightel] += jumpNorm;
2093  }//for i
2094 
2095  for(int64_t i = 0; i < nelements; i++){
2096  if (!AllCels[i]) continue;
2097  TPZCompElDisc * disc = dynamic_cast<TPZCompElDisc*>(AllCels[i]);
2098  if (!disc) continue;
2099  if(disc->Reference()->Dimension() != dim) continue;
2100  const STATE celJumpError = celJumps[i];
2101  //const STATE celJumpError = celJumps[i];
2102  if (abs(celJumpError) < eps){
2103  int64_t index;
2104  this->Discontinuous2Continuous(i, index);
2105  }//if
2106  }//for i
2107 
2108  this->CleanUpUnconnectedNodes();
2109  this->AdjustBoundaryElements();
2110 
2111 }//method
2112 
2113 void TPZCompMesh::AssembleError(TPZFMatrix<REAL> &estimator, int errorid){
2114  int64_t iel, i;
2115  const int64_t nelem = this->NElements();
2116  TPZManVector<REAL> locerror(7);
2117  TPZManVector<STATE> errorL(7), errorR(7);
2118 
2119  estimator.Resize(nelem, 7);
2120  estimator.Zero();
2121  for(iel=0; iel < nelem; iel++) {
2122  TPZCompEl *el = fElementVec[iel];
2123  if (!el) continue;
2124  TPZInterfaceElement * face = dynamic_cast<TPZInterfaceElement*>(el);
2125  if (face){
2126  errorL.Fill(0.); errorR.Fill(0.);
2127  face->ComputeErrorFace(errorid, errorL, errorR);
2128  int64_t n = errorL.NElements();
2129  if (errorR.NElements() > n) n = errorR.NElements();
2130  //if number of errors > 1 then resize matrix.
2131  //Method Resize keeps previous values and zero new values.
2132  //estimator.Resize(nelem, n);
2133  for(i = 0; i < errorL.NElements()-1; i++) {
2134  estimator(face->LeftElement()->Index(),i) += fabs(errorL[i]);
2135  }//for i
2136  for(i = 0; i < errorR.NElements()-1; i++){
2137  estimator(face->RightElement()->Index(),i) += fabs(errorR[i]);
2138  }//for i
2139  }//if
2140  else{
2141  locerror.Fill(0.);
2142  el->ComputeError(errorid, locerror);
2143  //if number of errors > 1 then resize matrix.
2144  //Method Resize keeps previous values and zero new values.
2145  //estimator.Resize(nelem, locerror.NElements());
2146  for(i = 0; i < locerror.NElements()-1; i++){
2147  estimator(iel, i) += locerror[i];
2148  }//for i
2149  // estimator.Print("Estimator", std::cout, EFormatted);
2150 
2151  }//else
2152  }
2153  for(iel=0; iel < nelem; iel++) {
2154  if(fabs(estimator(iel,5))>=0.0000001){ // Possivel erro ao plotar os indíces de efetividade
2155  estimator(iel,6) = estimator(iel,2)/estimator(iel,5) ;
2156  }
2157  }
2158 
2159 }
2160 
2162 TPZVec<STATE> TPZCompMesh::Integrate(const std::string &varname, const std::set<int> &matids)
2163 {
2164  // the postprocessed index of the varname for each material id
2165  std::map<int,int> variableids;
2166  int nvars = 0;
2167 
2168  std::map<int,TPZMaterial *>::iterator itmap;
2169  for (itmap = MaterialVec().begin(); itmap != MaterialVec().end(); itmap++) {
2170  if (matids.find(itmap->first) != matids.end()) {
2171  TPZMaterial *mat = itmap->second;
2172  TPZBndCond *bndcond = dynamic_cast<TPZBndCond *>(mat);
2173  int varindex = mat->VariableIndex(varname);
2174  if (varindex == -1 && bndcond) {
2175  mat = bndcond->Material();
2176  varindex= mat->VariableIndex(varname);
2177  }
2178  if (varindex == -1) {
2179  DebugStop();
2180  }
2181  variableids[itmap->first] = varindex;
2182  int nvarnew = mat->NSolutionVariables(variableids[itmap->first]);
2183  // the number of variables has to be the same for all materials
2184  if(nvars && nvars != nvarnew)
2185  {
2186  DebugStop();
2187  }
2188  nvars = nvarnew;
2189  }
2190  }
2191  TPZManVector<STATE,3> result(nvars,0.);
2192  int64_t nelem = NElements();
2193  for (int64_t el=0; el<nelem; el++) {
2194  TPZCompEl *cel = Element(el);
2195  if (!cel) {
2196  continue;
2197  }
2198  TPZGeoEl *gel = cel->Reference();
2199  if (!gel) {
2200  TPZManVector<STATE,3> locres;
2201  locres = cel->IntegrateSolution(varname, matids);
2202  if (locres.size() == nvars) {
2203  for (int iv = 0; iv<nvars; iv++) {
2204  result[iv] += locres[iv];
2205  }
2206  }
2207  else if(locres.size())
2208  {
2209  DebugStop();
2210  }
2211  }
2212  else
2213  {
2214  int matid = gel->MaterialId();
2215  if (matids.find(matid) == matids.end()) {
2216  continue;
2217  }
2218  TPZManVector<STATE,3> locres(nvars,0.);
2219  locres = cel->IntegrateSolution(variableids[matid]);
2220  if (locres.size() != nvars) {
2221  DebugStop();
2222  }
2223  for (int iv=0; iv<nvars; iv++)
2224  {
2225  result[iv] += locres[iv];
2226  }
2227  // std::cout << "el = " << el << " integral " << locres << " result " << result << std::endl;
2228  }
2229  }
2230  return result;
2231 }
2232 
2233 
2234 /*
2235 void TPZCompMesh::SaddlePermute()
2236 {
2237 
2238 #ifdef LOG4CXX
2239  {
2240  std::stringstream sout;
2241  sout<< "Implementando permutacao para problemas de ponto de sela"<< std::endl;
2242  LOGPZ_DEBUG(logger, sout.str().c_str());
2243  }
2244 #endif
2245  TPZVec<int64_t> permute;
2246  int numinternalconnects = NIndependentConnects();
2247  permute.Resize(numinternalconnects,0);
2248 
2249  TPZSubCompMesh *submesh = dynamic_cast<TPZSubCompMesh *> (this);
2250  if(submesh)
2251  {
2252  int nexternal = submesh->NConnects();
2253  numinternalconnects -= nexternal;
2254  }
2255  // else {
2256  // DebugStop();
2257  // }
2258 
2259  int jperm=0;
2260  int nel=ElementVec().NElements();
2261  for (int jel=0; jel<nel; jel++) {
2262 
2263  for (int64_t ip=0; ip<permute.NElements(); ip++) {
2264  permute[ip]=ip;
2265  }
2266 
2267  TPZCompEl *cel= ElementVec()[jel];
2268  // int idtroca=0;
2269  int eqmax=0;
2270  if(!cel)continue;
2271 // int ncon=elvec->NConnects();
2272  std::set<int> connects;
2273  cel->BuildConnectList(connects );
2274  // if(ncon==1) continue;
2275  int pressureconectindex = cel->PressureConnectIndex();
2276  if(pressureconectindex == -1) continue;
2277  int64_t eqpress=cel->Connect(pressureconectindex).SequenceNumber();
2278 
2279  for (std::set<int>::const_iterator it= connects.begin(); it != connects.end(); it++) {
2280 // for (int icon=0; icon< ncon-1; icon++) {
2281  if(*it == pressureconectindex) continue;
2282  TPZConnect &coel= fConnectVec[*it];
2283  if(coel.HasDependency()) continue;
2284  int eqflux=coel.SequenceNumber();
2285  eqmax = max(eqmax,eqflux);
2286  }
2287 
2288 
2289  if(eqpress<eqmax) {
2290 
2291  permute[eqpress]=eqmax;
2292 
2293  }
2294 
2295 
2296  for ( jperm = eqpress+1; jperm<=eqmax; jperm++) {
2297  permute[jperm]=jperm-1;
2298 
2299  }
2300 
2301 // #ifdef LOG4CXX
2302 // {
2303 // std::stringstream sout;
2304 // sout << "vetor SaddlePermute do elemento - "<<jel<< " - " <<permute;
2305 // LOGPZ_DEBUG(logger, sout.str().c_str());
2306 // }
2307 // #endif
2308 
2309  Permute(permute);
2310 
2311  }
2312 }
2313 */
2314 
2315 static void switchEq(int64_t eqsmall, int64_t eqlarge, TPZVec<int64_t> &permutegather, TPZVec<int64_t> &permutescatter)
2316 {
2317  int64_t eqkeep = permutegather[eqsmall];
2318  for (int64_t eq = eqsmall; eq< eqlarge; eq++) {
2319  permutegather[eq] = permutegather[eq+1];
2320  }
2321  permutegather[eqlarge] = eqkeep;
2322  for (int64_t eq = eqsmall; eq<= eqlarge; eq++) {
2323  permutescatter[permutegather[eq]] = eq;
2324  }
2325 }
2326 
2327 
2329 {
2330  TPZVec<int64_t> permutegather,permutescatter;
2331  int64_t numinternalconnects = NIndependentConnects();
2332  permutegather.Resize(numinternalconnects,0);
2333  permutescatter.Resize(numinternalconnects,0);
2334  for (int64_t i=0; i<numinternalconnects; i++) {
2335  permutegather[i] = i;
2336  permutescatter[i] = i;
2337  }
2338  int64_t numconnects = ConnectVec().NElements();
2339  int64_t numindepconnects = NIndependentConnects();
2340  if (numconnects==0) {
2341  return;
2342  }
2343  int minlagrange = 0;
2344  int maxlagrange = 0;
2345  for (int64_t ic=0; ic<numconnects; ic++) {
2346  TPZConnect &c = ConnectVec()[ic];
2347  if(c.HasDependency() || c.IsCondensed()) continue;
2348  if (c.SequenceNumber() < 0) {
2349  continue;
2350  }
2351  minlagrange = c.LagrangeMultiplier();
2352  maxlagrange = c.LagrangeMultiplier();
2353  break;
2354  }
2355  for (int ic=0; ic<numconnects; ic++) {
2356  TPZConnect &c = ConnectVec()[ic];
2357  if(c.HasDependency() || c.IsCondensed()) continue;
2358  int lagrange = c.LagrangeMultiplier();
2359  minlagrange = Min(lagrange, minlagrange);
2360  maxlagrange = Max(lagrange,maxlagrange);
2361  }
2362 
2363  int64_t nel = NElements();
2364  for (int lagr = minlagrange+1; lagr <= maxlagrange; lagr++)
2365  {
2366  for (int64_t el = 0; el<nel ; el++) {
2367  TPZCompEl *cel = ElementVec()[el];
2368  if (!cel) {
2369  continue;
2370  }
2371  int nc = cel->NConnects();
2372  if (nc <= 1) {
2373  continue;
2374  }
2375 
2376 #ifdef LOG4CXX
2377  if(logger->isDebugEnabled())
2378  {
2379  std::stringstream sout;
2380  sout << "el " << el << " Before renumbering : ";
2381  for (int ic=0; ic<nc; ic++) {
2382  TPZConnect &c = cel->Connect(ic);
2383  if (c.HasDependency() || c.IsCondensed()) {
2384  continue;
2385  }
2386  sout << permutescatter[c.SequenceNumber()] << "/" << (int)c.LagrangeMultiplier() << " ";
2387  }
2388  LOGPZ_DEBUG(logger, sout.str())
2389  }
2390 #endif
2391  // put all connects after the connect largest seqnum and lower lagrange number
2392  int64_t maxseq = -1;
2393  for (int ic=0; ic<nc ; ic++) {
2394  TPZConnect &c = cel->Connect(ic);
2395  if (c.HasDependency() || c.IsCondensed()) {
2396  continue;
2397  }
2398  int64_t eq = permutescatter[c.SequenceNumber()];
2399  if (!c.HasDependency() && c.LagrangeMultiplier() < lagr && eq > maxseq) {
2400  maxseq = eq;
2401  }
2402  }
2403  std::set<int64_t> seteq;
2404  for (int ic=0; ic<nc; ic++) {
2405  TPZConnect &c = cel->Connect(ic);
2406  if (c.HasDependency() || c.IsCondensed()) {
2407  continue;
2408  }
2409  int eq = permutescatter[c.SequenceNumber()];
2410  if (c.LagrangeMultiplier() == lagr && eq < maxseq) {
2411 #ifdef PZDEBUG
2412  if (c.HasDependency()) {
2413  DebugStop();
2414  }
2415 #endif
2416  seteq.insert(eq);
2417  }
2418  }
2419  std::set<int64_t>::reverse_iterator it;
2420  int64_t count = 0;
2421  for (it = seteq.rbegin(); it != seteq.rend(); it++) {
2422  int64_t eq = *it;
2423 #ifdef LOG4CXX
2424  if (logger->isDebugEnabled()) {
2425  std::stringstream sout;
2426  sout << "Switch ceq = " << eq << " with maxeq = " << maxseq-count;
2427  // sout << permute;
2428  LOGPZ_DEBUG(logger, sout.str())
2429  }
2430 #endif
2431  // we have to switch
2432  switchEq(eq, maxseq-count, permutegather, permutescatter);
2433  count++;
2434  }
2435 
2436 #ifdef PZDEBUG
2437  for (int64_t i=0; i<numinternalconnects; i++) {
2438  if (permutescatter[permutegather[i]] != i) {
2439  std::cout << "permutegather " << permutegather << std::endl;
2440  std::cout << "permutescatter " << permutescatter << std::endl;
2441  DebugStop();
2442  }
2443  }
2444  // put all connects after the connect largest seqnum and lower lagrange number
2445  maxseq = -1;
2446  for (int ic=0; ic<nc ; ic++) {
2447  TPZConnect &c = cel->Connect(ic);
2448  if (c.IsCondensed() || c.HasDependency()) {
2449  continue;
2450  }
2451  int64_t eq = permutescatter[c.SequenceNumber()];
2452  if (c.LagrangeMultiplier() < lagr && eq < maxseq) {
2453  maxseq = eq;
2454  }
2455  }
2456  for (int ic=0; ic<nc; ic++) {
2457  TPZConnect &c = cel->Connect(ic);
2458  if (c.HasDependency() || c.IsCondensed()) {
2459  continue;
2460  }
2461  int eq = permutescatter[c.SequenceNumber()];
2462  if (c.LagrangeMultiplier() == lagr && eq < maxseq) {
2463  // we have to switch
2464  DebugStop();
2465  }
2466  }
2467 #endif
2468 #ifdef LOG4CXX
2469  if(logger->isDebugEnabled())
2470  {
2471  std::stringstream sout;
2472  sout << "el " << el << " After renumbering : ";
2473  for (int ic=0; ic<nc; ic++) {
2474  TPZConnect &c = cel->Connect(ic);
2475  if (c.HasDependency() || c.IsCondensed()) {
2476  continue;
2477  }
2478  sout << permutescatter[c.SequenceNumber()] << "/" << (int)c.LagrangeMultiplier() << " ";
2479  }
2480  LOGPZ_DEBUG(logger, sout.str())
2481  }
2482 #endif
2483  }
2484 
2485  }
2486 #ifdef LOG4CXX2
2487  if (logger->isDebugEnabled())
2488  {
2489  std::stringstream sout;
2490  sout << "Saddle permute new permutation ";
2491  for (int i = 0; i< permutescatter.size(); i++) {
2492  int64_t jmax = i+10;
2493  if (jmax > permutescatter.size()) {
2494  jmax = permutescatter.size();
2495  }
2496  for (int64_t j=i; j< jmax; j++) {
2497  sout << permutescatter[j] << ' ';
2498  }
2499  sout << std::endl;
2500  i+= 10;
2501  }
2502  LOGPZ_DEBUG(logger, sout.str())
2503  }
2504 #endif
2505 
2506  Permute(permutescatter);
2507 #ifdef PZDEBUG
2508 
2509 #ifdef LOG4CXX
2510  if (logger->isDebugEnabled()) {
2511  LOGPZ_DEBUG(logger, "******************* AFTER PERMUTATION **************************")
2512  }
2513 #endif
2514 
2515  for (int64_t i=0L; i<numinternalconnects; i++) {
2516  permutegather[i] = i;
2517  permutescatter[i] = i;
2518  }
2519  for (int64_t el = 0L; el<nel ; el++) {
2520  TPZCompEl *cel = ElementVec()[el];
2521  if (!cel) {
2522  continue;
2523  }
2524  int nc = cel->NConnects();
2525  if (nc == 0) {
2526  continue;
2527  }
2528 #ifdef LOG4CXX
2529  if(logger->isDebugEnabled())
2530  {
2531  std::stringstream sout;
2532  sout << "el " << el << " Final numbering : ";
2533  for (int ic=0; ic<nc; ic++) {
2534  TPZConnect &c = cel->Connect(ic);
2535  if (c.HasDependency() || c.IsCondensed()) {
2536  continue;
2537  }
2538  sout << permutescatter[c.SequenceNumber()] << "/" << (int)c.LagrangeMultiplier() << " ";
2539  }
2540  LOGPZ_DEBUG(logger, sout.str())
2541  }
2542 #endif
2543  TPZConnect &c0 = cel->Connect(0);
2544  int minlagrange = c0.LagrangeMultiplier();
2545  int maxlagrange = c0.LagrangeMultiplier();
2546  for (int ic=0; ic<nc; ic++) {
2547  TPZConnect &c = cel->Connect(ic);
2548  if(c.HasDependency() || c.IsCondensed()) continue;
2549  int lagrange = c.LagrangeMultiplier();
2550  minlagrange = Min(lagrange, minlagrange);
2551  maxlagrange = Max(lagrange,maxlagrange);
2552  }
2553  for (int lagr = minlagrange+1; lagr <= maxlagrange; lagr++) {
2554  // put all connects after the connect largest seqnum and lower lagrange number
2555  int64_t maxseq = -1;
2556  for (int ic=0; ic<nc ; ic++) {
2557  TPZConnect &c = cel->Connect(ic);
2558  if (c.HasDependency() || c.IsCondensed()) {
2559  continue;
2560  }
2561  int64_t eq = permutescatter[c.SequenceNumber()];
2562  if (!c.HasDependency() && c.LagrangeMultiplier() < lagr && eq > maxseq) {
2563  maxseq = eq;
2564  }
2565  }
2566  for (int ic=0; ic<nc; ic++) {
2567  TPZConnect &c = cel->Connect(ic);
2568  if (c.HasDependency() || c.IsCondensed()) {
2569  continue;
2570  }
2571  int eq = permutescatter[c.SequenceNumber()];
2572  if (c.LagrangeMultiplier() == lagr && eq < maxseq) {
2573  // we have to switch
2574  cel->Print(std::cout);
2575  for (int iic=0; iic<nc; iic++) {
2576  cel->Connect(iic).Print(*this, std::cout);
2577  std::cout << "Destination seqnum " << permutegather[cel->Connect(iic).SequenceNumber()] << std::endl;
2578  }
2579  DebugStop();
2580  }
2581  }
2582  }
2583  }
2584 
2585 #endif
2586 }
2587 
2589 {
2590  TPZVec<int64_t> permute;
2591  int64_t numinternalconnects = NIndependentConnects();
2592  permute.Resize(numinternalconnects,0);
2593  for (int64_t i=0L; i<numinternalconnects; i++) {
2594  permute[i] = i;
2595  }
2596  int64_t nel = NElements();
2597  for (int64_t el = 0L; el<nel ; el++) {
2598  TPZCompEl *cel = ElementVec()[el];
2599  if (!cel) {
2600  continue;
2601  }
2602  unsigned char minlagrange = 255, maxlagrange = 0;
2603  int nc = cel->NConnects();
2604  for (int ic=0; ic<nc; ic++) {
2605  TPZConnect &c = cel->Connect(ic);
2606  int64_t seqnum = c.SequenceNumber();
2607  // if seqnum is larger than internal connects the equation is restrained, no permutation is necessary
2608  if (seqnum >= numinternalconnects) {
2609  continue;
2610  }
2611  unsigned char lagrange = c.LagrangeMultiplier();
2612  if (lagrange < minlagrange) {
2613  minlagrange = lagrange;
2614  }
2615  if (lagrange > maxlagrange) {
2616  maxlagrange = lagrange;
2617  }
2618  }
2619  for (unsigned char lagrange = minlagrange+1; lagrange <= maxlagrange; lagrange++) {
2620  int64_t maxeq = -1;
2621  for (int ic=0; ic<nc ; ic++) {
2622  TPZConnect &c = cel->Connect(ic);
2623  if (c.SequenceNumber() >= numinternalconnects) {
2624  continue;
2625  }
2626  if (c.LagrangeMultiplier() < lagrange) {
2627  int64_t origeq = c.SequenceNumber();
2628  if(maxeq < permute[origeq])
2629  {
2630  maxeq = permute[origeq];
2631  }
2632  }
2633  }
2634  if (maxeq < 0) {
2635  continue;
2636  }
2637  std::set<int64_t> lagrangeseqnum;
2638  for (int ic=nc-1; ic>=0 ; ic--) {
2639  TPZConnect &c = cel->Connect(ic);
2640  int clagrange = c.LagrangeMultiplier();
2641  int64_t ceqnum = c.SequenceNumber();
2642  if (ceqnum >= numinternalconnects) {
2643  continue;
2644  }
2645  int64_t ceq = permute[ceqnum];
2646  if (clagrange == lagrange && ceq < maxeq) {
2647  lagrangeseqnum.insert(ceqnum);
2648  }
2649  }
2650  std::set<int64_t>::reverse_iterator it;
2651  int64_t count = 0;
2652  for (it = lagrangeseqnum.rbegin(); it != lagrangeseqnum.rend(); it++) {
2653  int64_t ceq = permute[*it];
2654  ModifyPermute(permute, ceq, maxeq-count);
2655 #ifdef LOG4CXX
2656  if(logger->isDebugEnabled())
2657  {
2658  std::stringstream sout;
2659  sout << "Switch ceq = " << ceq << " with maxeq = " << maxeq-count;
2660 // sout << permute;
2661  LOGPZ_DEBUG(logger, sout.str())
2662  }
2663 #endif
2664  count++;
2665  }
2666 #ifdef LOG4CXX
2667  if(logger->isDebugEnabled())
2668  {
2669  std::stringstream sout;
2670  sout << "Resequence for element " << el << std::endl;
2671  for (int ic=0; ic<nc; ic++) {
2672  TPZConnect &c = cel->Connect(ic);
2673  c.Print(*this,sout);
2674  if (c.SequenceNumber() < numinternalconnects) {
2675  sout << "New seqnum = " << permute[c.SequenceNumber()] << std::endl;
2676  }
2677  else
2678  {
2679  sout << "Connect with restraint, seqnum " << c.SequenceNumber() << std::endl;
2680  }
2681  }
2682  LOGPZ_DEBUG(logger, sout.str())
2683  }
2684 #endif
2685  }
2686  }
2687 #ifdef LOG4CXX2
2688  if(logger->isDebugEnabled())
2689  {
2690  for (int64_t el=0L; el<nel; el++) {
2691  TPZCompEl *cel = ElementVec()[el];
2692  if (!cel) {
2693  continue;
2694  }
2695  int nc = cel->NConnects();
2696  std::stringstream sout;
2697  sout << "Resequence for element " << el << std::endl;
2698  for (int ic=0; ic<nc; ic++) {
2699  TPZConnect &c = cel->Connect(ic);
2700  c.Print(*this,sout);
2701  if (c.SequenceNumber() < numinternalconnects) {
2702  sout << "New seqnum = " << permute[c.SequenceNumber()] << std::endl;
2703  }
2704  }
2705  LOGPZ_DEBUG(logger, sout.str())
2706  }
2707  }
2708 #endif
2709 #ifdef LOG4CXX
2710  if (logger->isDebugEnabled())
2711  {
2712  std::stringstream sout;
2713  sout << "Saddle permute old permutation ";
2714  for (int i = 0; i< permute.size(); i++) {
2715  int64_t jmax = i+10;
2716  if (jmax > permute.size()) {
2717  jmax = permute.size();
2718  }
2719  for (int64_t j=i; j< jmax; j++) {
2720  sout << permute[j] << ' ';
2721  }
2722  sout << std::endl;
2723  i+= 10;
2724  }
2725  LOGPZ_DEBUG(logger, sout.str())
2726  }
2727 #endif
2728  Permute(permute );
2729 
2730 }
2731 
2733 void TPZCompMesh::ModifyPermute(TPZVec<int64_t> &permute, int64_t lagrangeq, int64_t maxeq)
2734 {
2735  int64_t neq = permute.size();
2736 #ifdef PZDEBUG
2737  if (lagrangeq < 0 || lagrangeq >= neq || maxeq < 0 || maxeq >= neq) {
2738  DebugStop();
2739  }
2740 #endif
2741  // find the equation which maps to lagrangeq
2742  //int lagrangeqindex = permuteinv[lagrangeq];
2743  TPZVec<int64_t> accpermute(neq,0),input(permute);
2744  for (int64_t i=0; i<neq; i++) {
2745  accpermute[i] = i;
2746  }
2747 
2748  int64_t lagrangeqindex = lagrangeq;
2749 
2750  // this equation should never be sent forwards
2751  if (accpermute[lagrangeqindex] > lagrangeq) {
2752  DebugStop();
2753  }
2754 
2755  accpermute[lagrangeqindex] = maxeq;
2756  int64_t index = lagrangeqindex+1;
2757  while (index < neq && (accpermute[index] <= maxeq || accpermute[index] < index)) {
2758  accpermute[index] = accpermute[index]-1;
2759  index++;
2760  }
2761  for (int64_t i=0; i<neq; i++) {
2762  permute[i] = accpermute[input[i]];
2763  }
2764 
2765 #ifdef PZDEBUG22
2766  {
2767  std::set<int64_t> acc;
2768  for (int64_t i=0; i<neq; i++) {
2769  acc.insert(permute[i]);
2770  }
2771  if (acc.size() != neq) {
2772  std::cout << "input " << input << std::endl;
2773  std::cout << "accpermute " << accpermute << std::endl;
2774  std::cout << "permute " << permute << std::endl;
2775  DebugStop();
2776  }
2777  }
2778 #endif
2779 }
2780 
2782 void TPZCompMesh::BuildCornerConnectList(std::set<int64_t> &connectindexes) const
2783 {
2784  int64_t nel = NElements();
2785  for (int64_t el=0; el<nel ; el++) {
2786  TPZCompEl *cel = ElementVec()[el];
2787  if (!cel) {
2788  continue;
2789  }
2790  cel->BuildCornerConnectList(connectindexes);
2791  }
2792 }
2793 
2795 
2796  TPZStack<TPZCompMesh *> s1, s2;
2797  int64_t pos1=0, pos2, comind;
2798  TPZCompMesh *father = FatherMesh();
2799  s1.Push(this);
2800  while (father){
2801  s1.Push((father));
2802  pos1++;
2803  father = s1[pos1]->FatherMesh();
2804  }
2805  pos2 = 0;
2806  s2.Push(mesh);
2807  father = mesh->FatherMesh();
2808  while (father){
2809  s2.Push(father);
2810  pos2++;
2811  father = s2[pos2]->FatherMesh();
2812  }
2813  if (s1[pos1] != s2[pos2]) return 0;
2814  comind=0; //The first mesh is common for all submeshes
2815  for (; pos1>=0 && pos2>=0; pos1--, pos2--) {
2816  if((s1[pos1])!=(s2[pos2])) {
2817  comind=pos1+1;
2818  return (s1[comind]);
2819  }
2820  }
2821  return (pos1 >=0 ) ? (s1[pos1+1]) : s2[pos2+1];
2822 }
2823 
2827 {
2828  if(fSolN.Rows() != fSolution.Rows() || fSolN.Cols() != fSolution.Cols())
2829  {
2831  }
2832  fSolN += mult*fSolution;
2833  fSolution = fSolN;
2834  int64_t nel = NElements();
2835  for(int64_t el = 0; el<nel; el++)
2836  {
2837  TPZCompEl *cel = Element(el);
2838  TPZSubCompMesh *subc = dynamic_cast<TPZSubCompMesh *>(cel);
2839  if(subc) subc->UpdatePreviousState(mult);
2840  }
2841 }
2842 
2843 
2844 template class TPZRestoreClass<TPZCompMesh>;
2845 
int64_t NElements() const
Number of computational elements allocated.
Definition: pzcmesh.h:169
virtual TPZVec< STATE > IntegrateSolution(int var) const
Compute the integral of a variable.
Definition: pzcompel.cpp:1054
TPZGeoMesh * Reference() const
Returns a pointer to the geometrical mesh associated.
Definition: pzcmesh.h:209
void SetAllCreateFunctionsDiscontinuous()
Create discontinuous approximation spaces.
Contains TPZShapeTetra class which implements the shape functions of a tetrahedral element...
This class performs a series of consistency tests on geometric transformations between elements...
Definition: pzcheckgeom.h:16
int fDimModel
Definition: pzcmesh.h:86
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
void SetElementSolution(int64_t i, TPZVec< STATE > &sol)
Set a ith element solution, expanding the element-solution matrix if necessary.
Definition: pzcmesh.cpp:1533
virtual TPZCompMesh * FatherMesh() const
Get the father meshes stack.
Definition: pzcmesh.h:337
void AdjustBoundaryElements()
Will refine the elements associated with a boundary condition till there are no elements constrained ...
Definition: pzcmesh.cpp:1423
void CreateInterfaces(bool BetweenContinuous=false)
Create interfaces between this and its neighbours.
virtual void BuildCornerConnectList(std::set< int64_t > &connectindexes) const =0
adds the connect indexes associated with base shape functions to the set
TPZCompMesh * Reference() const
Returns the currently loaded computational grid.
Definition: pzgmesh.h:167
int Position(const int block_diagonal) const
Returns the position of first element block dependent on matrix diagonal.
Definition: pzblock.h:177
void IncrementElConnected()
Increment fNElConnected.
Definition: pzconnect.h:260
void SetElementGraph(TPZVec< int64_t > &elgraph, TPZVec< int64_t > &elgraphindex)
This method declares the element graph to the object.
Implements computational element and a side. Computational Element.
Definition: pzcompel.h:632
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
Represents a set of shape functions associated with a computational element/side. Computational Eleme...
Definition: pzconnect.h:30
TPZGeoNode * NodePtr(int i) const
Returns a pointer to the ith node of the element.
Definition: pzgeoel.cpp:2566
void SetName(const std::string &nm)
Set the mesh name.
Definition: pzcmesh.cpp:232
Contains declaration of TPZGeoElSide class which represents an element and its side, and TPZGeoElSideIndex class which represents an TPZGeoElSide index.
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
virtual MElementType Type()
Return the type of the element.
Definition: pzcompel.cpp:195
Contains declaration of TPZGeoElRefLess class which implements the mapping between the master element...
TPZFMatrix< STATE > fSolution
Solution vector.
Definition: pzcmesh.h:72
int MaterialId() const
Returns the material index of the element.
Definition: pzgeoel.h:250
void Unwrap()
put the elements in the element group back in the mesh and delete the element group ...
virtual void BuildConnectList(std::set< int64_t > &indepconnectlist, std::set< int64_t > &depconnectlist)
Builds the list of all connectivities related to the element including the connects pointed to by dep...
Definition: pzcompel.cpp:435
TPZAdmChunkVector< TPZCompEl * > fElementVec
List of pointers to elements.
Definition: pzcmesh.h:60
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
Implements renumbering for elements of a mesh. Utility.
Definition: pzmetis.h:17
void ConvertGraph(TPZVec< int64_t > &elgraph, TPZVec< int64_t > &elgraphindex, TPZManVector< int64_t > &nodegraph, TPZManVector< int64_t > &nodegraphindex)
Will convert an element graph defined by elgraph and elgraphindex into a node graph defined by nodegr...
Contains declaration of TPZGeoNode class which defines a geometrical node.
REAL Coord(int i) const
Returns i-th coordinate of the current node.
Definition: pzgnode.h:99
int64_t NEquations()
This computes the number of equations associated with non-restrained nodes.
Definition: pzcmesh.cpp:721
void ConvertDiscontinuous2Continuous(REAL eps, int opt, int dim, TPZVec< STATE > &celJumps)
This method compute the jump solution of interface and convert discontinuous elements with jump less ...
Definition: pzcmesh.cpp:2052
int64_t NIndependentConnects()
Number of independent connect objects.
Definition: pzcmesh.cpp:1080
void ProjectSolution(TPZFMatrix< STATE > &projectsol)
void SaddlePermute()
Put the sequence number of the pressure connects after the seq number of the flux connects...
Definition: pzcmesh.cpp:2328
int64_t fNmeshes
Definition: pzcmesh.h:98
static REAL Distance(TPZVec< REAL > &centel, TPZVec< REAL > &centface)
Definition: pzgeoel.cpp:936
Contains declaration of TPZCompEl class which defines the interface of a computational element...
Templated vector implementation.
void WritePointers(const TPZVec< TPZAutoPointer< T >> &vec)
Definition: TPZStream.h:214
Contains the declaration of the TPZElementGroup class, which implements an computational element whic...
int ClassId() const override
Returns the unique identifier for reading/writing objects to streams.
Definition: pzcmesh.cpp:1969
TPZGeoElSide Father2() const
returns the father/side pair which contains the set of points associated with this/side ...
virtual TPZTransform< REAL > BuildTransform2(int side, TPZGeoEl *father, TPZTransform< REAL > &t)
Returns the transformation which maps the parameter side of the element/side into the parameter spac...
Definition: pzgeoel.cpp:386
virtual int Dimension() const =0
Dimension of the element.
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
void SetBlocks(TPZBlock< TVar > &row, TPZBlock< TVar > &col, int nvar, int nrowblocks, int ncolblocks)
This operation will reset the matrix to zero with no rows defined.
Definition: pztransfer.cpp:85
Defines PZError.
virtual int NCornerNodes() const =0
Returns the number of corner nodes of the element.
TPZVec< STATE > Integrate(const std::string &varname, const std::set< int > &matids)
Integrate the variable name over the mesh.
Definition: pzcmesh.cpp:2162
Contains the TPZRefQuad class which implements the uniform refinement of a geometric quadrilateral el...
REAL LesserEdgeOfMesh()
Definition: pzcmesh.cpp:1849
Contains declaration of TPZInterpolationSpace class which implements the interface for interpolated c...
virtual void RemoveSideRestraintsII(MInsertMode mode)
Delete the restraints on the nodes of the connected elements if necessary.
Definition: pzintel.cpp:1286
virtual int Dimension() const =0
Returns the integrable dimension of the material.
virtual int NStateVariables() const =0
Returns the number of state variables associated with the material.
Contains the declaration of the TPZMultiphysicsElement class. This class is abstract.
Declarates the TPZBlock<REAL>class which implements block matrices.
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
Class which groups elements to characterize dense matrices.
static int fatherside[8][21]
Definition: pzrefprism.cpp:303
size_t NMaterials() const
Number of materials.
Definition: pzcmesh.h:172
Defines enum MElementType and contains the implementation of MElementType_NNodes(...) functions.
TinyFad< 8, T > abs(const TinyFad< 8, T > &in)
Definition: tinyfadeight.h:846
clarg::argString input("-if", "input file", "cube1.txt")
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
int64_t NElements() const
Access method to query the number of elements of the vector.
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
Implements a chunk vector with free store administration. Utility.
Definition: TPZStream.h:39
void BuildTransferMatrix(TPZCompElDisc &coarsel, TPZTransfer< STATE > &transfer)
std::string & Name()
Definition: pzgmesh.h:135
void ComputeFillIn(int64_t resolution, TPZFMatrix< REAL > &fillin)
This method will fill the matrix passed as parameter with a representation of the fillin of the globa...
Definition: pzcmesh.cpp:1869
virtual int NSides() const =0
Returns the number of connectivities of the element.
static TPZSavable * GetInstance(const int64_t &objId)
void ReadPointers(TPZVec< TPZAutoPointer< T >> &vec)
Definition: TPZStream.h:305
virtual int SideDimension(int side) const =0
Return the dimension of side.
virtual void BuildCornerConnectList(std::set< int64_t > &connectindexes) const
adds the connect indexes associated with base shape functions to the set
Definition: pzcmesh.cpp:2782
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzblock.cpp:504
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
void CleanUp()
Delete all the dynamically allocated data structures.
Definition: pzcmesh.cpp:161
Contains the TPZTransfer class which implements a rectangular sparse block matrix.
void LoadReferences()
Map this grid in the geometric grid.
Definition: pzcmesh.cpp:482
virtual void LoadElementReference()
Loads the geometric element reference.
Definition: pzcompel.cpp:601
void RemakeAllInterfaceElements()
Deletes all interfaces and rebuild them all.
Definition: pzcmesh.cpp:1281
void ComputeElGraph(TPZStack< int64_t > &elgraph, TPZVec< int64_t > &elgraphindex)
Computes the connectivity graph of the elements, as appropriate for the TPZRenumbering class...
Definition: pzcmesh.cpp:1091
REAL DeltaX()
Definition: pzcmesh.cpp:1810
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
Contains the TPZRefPyramid class which implements the uniform refinement of a geometric hexahedral el...
virtual void TransferMultiphysicsElementSolution()
Definition: pzcompel.h:415
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Definition: pzlog.h:95
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
void error(char *string)
Definition: testShape.cc:7
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
int Dim() const
Returns matrix dimension pointed by block.
Definition: pzblock.h:180
void BuildTransferMatrix(TPZInterpolationSpace &coarsel, TPZTransform<> &t, TPZTransfer< STATE > &transfer)
Accumulates the transfer coefficients between the current element and the coarse element into the t...
virtual std::string TypeName() const
Returns the type of the element as a string.
Definition: pzgeoel.h:262
Contains TPZShapeCube class which implements the shape functions of a hexaedral element.
void GetElementPatch(TPZVec< int64_t > nodtoelgraph, TPZVec< int64_t > nodtoelgraphindex, TPZStack< int64_t > &elgraph, TPZVec< int64_t > &elgraphindex, int64_t elind, TPZStack< int64_t > &patch)
Gives the element patch.
Definition: pzcmesh.cpp:1599
void BuildTransferMatrix(TPZCompMesh &coarsemesh, TPZTransfer< STATE > &transfer)
Builds the transfer matrix from the current grid to the coarse grid.
Definition: pzcmesh.cpp:883
static TPZAutoPointer< TPZSavable > GetAutoPointer(const int64_t &objId)
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
virtual TPZMaterial * Material() const
Identify the material object associated with the element.
Definition: pzcompel.cpp:959
unsigned char LagrangeMultiplier() const
Access method to return the indication whether the connect is associated with a pressure lagrange mul...
Definition: pzconnect.h:230
void HigherLevelElementList(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated, int removeduplicates)
Returns all connected elements which have level higher to the current element.
Definition: pzcompel.cpp:761
TPZAutoPointer< TPZGeoMesh > fGMesh
Autopointer to the geometric mesh used in case the user has passed an autopointer.
Definition: pzcmesh.h:54
virtual void AutoBuildContDisc(const TPZVec< TPZGeoEl *> &continuous, const TPZVec< TPZGeoEl *> &discontinuous)
Creates the computational elements, and the degree of freedom nodes.
Definition: pzcmesh.cpp:326
TPZCompEl * CreateCompEl(TPZGeoEl *gel, int64_t &index)
Create a computational element based on the geometric element.
Definition: pzcmesh.h:462
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
Contains the TPZRefPrism class which implements the uniform refinement of a geometric prism element...
TPZCompEl * Element(int64_t iel)
Definition: pzcmesh.h:185
void LoadSolution(const TPZFMatrix< STATE > &sol)
Given the solution of the global system of equations, computes and stores the solution for the restri...
Definition: pzcmesh.cpp:441
void Resize(const int newsize)
Increase the size of the chunk vector.
Definition: pzadmchunk.h:280
void Unwrap()
unwrap the condensed element from the computational element and delete the condensed element ...
virtual TPZCompEl * Clone(TPZCompMesh &mesh) const =0
Method for creating a copy of the element.
int fDefaultOrder
Default order for all elements of this mesh.
Definition: pzcmesh.h:89
Contains declaration of TPZMesh class which defines a geometrical mesh and contains a corresponding l...
int64_t SequenceNumber() const
Returns the Sequence number of the connect object.
Definition: pzconnect.h:158
void GetRefPatches(std::set< TPZGeoEl *> &grpatch)
Gives all patches of the mesh.
Definition: pzcmesh.cpp:1558
virtual void ComputeErrorFace(int errorid, TPZVec< STATE > &errorL, TPZVec< STATE > &errorR)
ComputeError computes the element error estimator.
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
void RemoveInterfaces()
Remove interfaces connected to this element.
virtual void Skyline(TPZVec< int64_t > &skyline)
This method computes the skyline of the system of equations.
Definition: pzcmesh.cpp:787
void Reset()
Reset the data of the connect.
Definition: pzconnect.h:113
void EqualLevelElementList(TPZStack< TPZCompElSide > &elsidevec, int onlyinterpolated, int removeduplicates)
Returns all connected elements which have equal level to the current element This method will not put...
Definition: pzcompel.cpp:771
virtual void Divide(int64_t index, TPZVec< int64_t > &subindex, int interpolate=0)
Divide the computational element. If interpolate = 1, the solution is interpolated to the sub element...
Definition: pzcompel.cpp:410
virtual void LoadSolution()
Loads the solution within the internal data structure of the element.
Definition: pzcompel.cpp:200
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
Contains declaration of.
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
const T & Max(const T &a, const T &b)
Returns the maximum value between a and b.
Definition: pzreal.h:724
int64_t Index() const
Returns the index of the element within the element vector of the mesh.
Definition: pzgeoel.h:730
virtual bool HasMaterial(const std::set< int > &materialids) const
Verifies if the material associated with the element is contained in the set.
Definition: pzcompel.cpp:968
Contains the TPZGeoTetrahedra class which implements the geometry of a tetrahedral element...
TPZCompMesh * Clone() const
Clone this mesh.
Definition: pzcmesh.cpp:1720
Implements a group of computational elements as a mesh and an element. Computational Mesh...
Definition: pzsubcmesh.h:36
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzcmesh.cpp:2006
Contains declaration of TPZInterfaceElement class which computes the contribution over an interface b...
virtual void ExpandSolution()
Adapt the solution vector to new block dimensions.
Definition: pzcmesh.cpp:396
void SetSequenceNumber(int64_t i)
Set the sequence number for the global system of equations of the connect object. ...
Definition: pzconnect.h:165
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZCompMesh * Mesh() const
Return a pointer to the grid of the element.
Definition: pzcompel.cpp:288
Contains declaration of TPZConnect class which represents a set of shape functions associated with a ...
Free store vector implementation.
int Dimension() const
Returns the dimension of the simulation.
Definition: pzcmesh.h:148
void Coarsen(TPZVec< int64_t > &elements, int64_t &index, bool CreateDiscontinuous=false)
Create a computational element father for the comp. elem. into elements.
Definition: pzcmesh.cpp:1164
virtual void ComputeNodElCon()
Compute the number of elements connected to each connect object.
Definition: pzcmesh.cpp:668
int NumInterfaces()
Returns number of TPZInterfaceElement pointing to this.
Definition: pzgeoel.h:94
void Divide(int64_t index, TPZVec< int64_t > &subindex, int interpolate=0)
Divide the element corresponding to index.
Definition: pzcmesh.cpp:1146
void SaddlePermute2()
Definition: pzcmesh.cpp:2588
TPZMaterial * FindMaterial(int id)
Find the material with identity id.
Definition: pzcmesh.cpp:297
virtual int64_t PutinSuperMesh(int64_t local, TPZCompMesh *super)
Put an local connection in the supermesh - Supermesh is one mesh who contains the analised submesh...
Definition: pzcmesh.cpp:1512
TPZCompEl * Reference() const
Return a pointer to the element referenced by the geometric element.
Definition: pzgeoel.h:750
unsigned int NShape() const
Definition: pzconnect.h:151
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
int HasDependency() const
Returns whether exist dependecy information.
Definition: pzconnect.h:292
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
void GetNodeToElGraph(TPZVec< int64_t > &nodtoelgraph, TPZVec< int64_t > &nodtoelgraphinde, TPZStack< int64_t > &elgraph, TPZVec< int64_t > &elgraphindexx)
Gives the conects graphs.
Definition: pzcmesh.cpp:1576
TPZFMatrix< STATE > fElementSolution
Solution vectors organized by element.
Definition: pzcmesh.h:83
virtual void ComputeError(int errorid, TPZVec< REAL > &error)
ComputeError computes the element error estimator.
Definition: pzcompel.h:353
void SetAllCreateFunctionsContinuous()
Create continuous approximation spaces.
const TPZBlock< STATE > & Block() const
Access the block structure of the solution vector.
Definition: pzcmesh.h:213
void Read(TPZStream &buf, void *context) override
read objects from the stream
static int GetgOrder()
Set the default value of the interpolation order.
Definition: pzcompel.h:830
int64_t NConnects() const
Number of connects allocated including free nodes.
Definition: pzcmesh.h:163
TPZCompMesh(TPZGeoMesh *gr=0)
Constructor from geometrical mesh.
Definition: pzcmesh.cpp:63
Contains declaration of TPZCompelDisc class which implements a computational element for discontinuou...
TPZGeoMesh * fReference
Geometric grid to which this grid refers.
Definition: pzcmesh.h:51
TPZAdmChunkVector< TPZConnect > & ConnectVec()
Return a reference to the connect pointers vector.
Definition: pzcmesh.h:198
TPZFMatrix< STATE > fSolN
Solution at previous state.
Definition: pzcmesh.h:76
Class which implements an element which condenses the internal connects.
int64_t Index() const
Returns element index of the mesh fELementVec list.
Definition: pzcompel.h:821
Contains the TPZMaterialData class which implements an interface between TPZCompEl::CalcStiff and TPZ...
const T & Min(const T &a, const T &b)
Returns the minimum value between a and b.
Definition: pzreal.h:729
Contains TPZShapePoint class which implements the shape function associated with a point...
virtual int PreferredSideOrder(int iside)=0
Returns the preferred order of the polynomial along side iside.
int SetNBlocks(const int num_of_blocks)
Sets number of blocks on diagonal matrix.
Definition: pzblock.cpp:91
Contains declaration of TPZCompMesh class which is a repository for computational elements...
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZBlock< STATE > fBlock
Block structure to right construction of the stiffness matrix and load vector.
Definition: pzcmesh.h:80
virtual TPZConnect & Connect(int i) const
Returns a pointer to the ith node.
Definition: pzcompel.cpp:298
void CompactDataStructure(CompactScheme type=CompactScheme::ALWAYS)
Sets the method to compact the data structure based on the.
Definition: pzadmchunk.h:222
int BandWidth()
This method computes the bandwidth of the system of equations.
Definition: pzcmesh.cpp:747
Contains the TPZGeoCube class which implements the geometry of hexahedra element. ...
void InterpolateSolution(TPZInterpolationSpace &coarsel)
Interpolates the solution into the degrees of freedom nodes from the degrees of freedom nodes from th...
void SetDimModel(int dim)
Set de dimension of the domain of the problem.
Definition: pzcmesh.h:139
unsigned char NState() const
Number of state variables associated with the connect.
Definition: pzconnect.h:146
void UpdatePreviousState(STATE mult=1.0)
Definition: pzcmesh.cpp:2826
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
Contains the TPZGeoPoint class which implements the geometry of a point element or 0-D element...
std::map< int,TPZMaterial *> & MaterialVec()
Returns a reference to the material pointers vector.
Definition: pzcmesh.h:203
Computes the contribution over an interface between two discontinuous elements. Computational Element...
Contains TPZMatrix<TVar>class, root matrix class.
virtual int HasSubElement() const =0
Return 1 if the element has subelements.
Contains TPZShapePrism class which implements the shape functions of a prism element.
Implements the interfaces for TPZCompElDisc, TPZInterfaceElement and TPZInterpolatedElement. Computational element.
int NBlocks() const
Returns number of blocks on diagonal.
Definition: pzblock.h:165
TPZBlock< STATE > fSolutionBlock
Block structure of the solution vector ????
Definition: pzcmesh.h:68
virtual int NConnects() const =0
Returns the number of nodes of the element.
int NeighbourExists(int side, const TPZGeoElSide &gel)
Returns 1 if gel is a neighbour of the element along side.
Definition: pzgeoel.cpp:226
Contains the TPZRefTriangle class which implements the uniform refinement of a geometric triangular e...
Contains declaration of TPZSubCompMesh class which implements a group of computational elements as a ...
static void switchEq(int64_t eqsmall, int64_t eqlarge, TPZVec< int64_t > &permutegather, TPZVec< int64_t > &permutescatter)
Definition: pzcmesh.cpp:2315
int InsertMaterialObject(TPZMaterial *mat)
Insert a material object in the datastructure.
Definition: pzcmesh.cpp:287
REAL LesserEdgeOfEl()
Will return the smallest distance between two nodes of the reference element.
Definition: pzcompel.cpp:697
Contains the TPZRefTetrahedra class which implements the uniform refinement of a geometric tetrahedra...
virtual void CleanUpUnconnectedNodes()
Delete the nodes which have no elements connected to them.
Definition: pzcmesh.cpp:498
virtual int Dimension() const =0
Returns the dimension of the element.
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
void AssembleError(TPZFMatrix< REAL > &estimator, int errorid)
Assemble the vector with errors estimators.
Definition: pzcmesh.cpp:2113
virtual ~TPZCompMesh()
Simple Destructor.
Definition: pzcmesh.cpp:134
int Dimension()
Get Dimension.
Definition: pzgmesh.h:291
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: pzfmatrix.cpp:2205
TPZGeoEl * Reference() const
Return a pointer to the corresponding geometric element if such exists, return 0 otherwise.
Definition: pzcompel.cpp:1137
Implements a geometric node in the pz environment. Geometry.
Definition: pzgnode.h:31
TPZCompEl * LeftElement() const
Returns the left element from the element interface.
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
TPZCompEl * SubElement(int64_t sub) const
Returns the "sub" subelement.
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
virtual REAL CompareMesh(int var, char *matname)
This method will initiate the comparison between the current computational mesh and the mesh which is...
Definition: pzcmesh.cpp:1522
void ResetElConnected()
Initialize with zero fNElConnected.
Definition: pzconnect.h:258
void ResetReference()
Reset the element referenced by the geometric element to NULL.
Definition: pzgeoel.h:431
Implements computational mesh. Computational Mesh.
Definition: pzcmesh.h:47
int Size(const int block_diagonal) const
Returns block dimension.
Definition: pzblock.h:171
TPZAdmChunkVector< TPZCompEl * > & ElementVec()
Returns a reference to the element pointers vector.
Definition: pzcmesh.h:183
void BuildTransferMatrixDesc(TPZCompMesh &transfermesh, TPZTransfer< STATE > &transfer)
To discontinuous elements.
Definition: pzcmesh.cpp:941
virtual void Print(std::ostream &out=std::cout) const
Prints element data.
Definition: pzcompel.cpp:322
void PRefine(int order) override
Changes the interpolation order of a side. Updates all constraints and block sizes ...
Definition: pzintel.cpp:1592
int NElConnected() const
Returns fNElConnected.
Definition: pzconnect.h:264
void BuildMesh(TPZCompMesh &cmesh, const std::set< int > &MaterialIDs) const
Creates the computational elements, and the degree of freedom nodes.
void SetMatrix(TPZMatrix< TVar > *const other)
Changes pointer to other.
Definition: pzblock.h:51
Contains declaration of TPZInterpolatedElement class which implements computational element of the in...
Contains the TPZRefLinear class which implements the uniform refinement of a geometric linear element...
void SetReference(TPZCompMesh *ref)
Sets the reference of the geometric grid to ref.
Definition: pzgmesh.h:161
int Id() const
Definition: TPZMaterial.h:170
virtual TPZCompMesh * CommonMesh(TPZCompMesh *mesh)
Gives the commom father mesh of the specified mesh and the current submesh.
Definition: pzcmesh.cpp:2794
void EvaluateError(std::function< void(const TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv)> fp, bool store_error, TPZVec< REAL > &errorSum)
Evaluates the error given the two vectors of the analised parameters.
Definition: pzcmesh.cpp:1395
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
TPZCompMesh & operator=(const TPZCompMesh &copy)
copy the content of the mesh
Definition: pzcmesh.cpp:1683
TPZCreateApproximationSpace fCreate
The object which defines the type of space being created.
Definition: pzcmesh.h:92
int64_t NIndexes() const
Returns the number of clustered subelements.
int Exists() const
Definition: pzgeoelside.h:175
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
virtual void SetDegree(int degree)
Assigns the degree of the element.
std::string fName
Grid name for model identification.
Definition: pzcmesh.h:57
void SetFree(int index)
Indicate an element as free.
Definition: pzadmchunk.h:199
pthread_cond_t cond
Definition: numatst.cpp:318
Contains TPZShapePiram class which implements the shape functions of a pyramid element.
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
Contains declaration of the abstract TPZStream class. TPZStream defines the interface for saving and ...
Contains the TPZGeoPyramid class which implements the geometry of pyramid element.
void CopyMaterials(TPZCompMesh &mesh) const
Copies the materials of this mesh to the given mesh.
Definition: pzcmesh.cpp:1797
Contains the TPZTransform<> class which implements an affine transformation between points in paramet...
std::map< int, TPZMaterial *> fMaterialVec
Map of pointers to materials.
Definition: pzcmesh.h:65
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZCompEl * RightElement() const
Returns the right element from the element interface.
void ModifyPermute(TPZVec< int64_t > &permute, int64_t lagrangeq, int64_t maxeq)
Modify the permute vector swapping the lagrangeq with maxeq and shifting the intermediate equations...
Definition: pzcmesh.cpp:2733
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
Contains the TPZGeoPrism class which implements the geometry of a prism element.
clarg::argInt porder("-porder", "polinomial order", 1)
Contains declaration of TPZGeoElement class which implements a generic geometric element with a unifo...
This class implements a discontinuous element (for use with discontinuous Galerkin). Computational Element.
Definition: TPZCompElDisc.h:35
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
void InitializeBlock()
Resequence the block object, remove unconnected connect objects and reset the dimension of the soluti...
Definition: pzcmesh.cpp:391
Implements rectangular matrix which extends a solution vector of the coarse mesh to a solution vector...
Definition: pzcmesh.h:33
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: pzadmchunk.h:118
REAL MaximumRadiusOfEl()
Will return the maximum distance between the nodes of the reference element.
Definition: pzcompel.cpp:671
void ProjectSolution(TPZFMatrix< STATE > &projectsol)
Definition: pzcmesh.cpp:1923
virtual void Print(std::ostream &out=std::cout)
Print all relevant data of the element to cout.
Definition: pzgeoel.cpp:238
void EvaluateInterfaceJump(TPZSolVec &jump, int opt)
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzfmatrix.cpp:2225
void ResetReference()
Resets all load references in elements and nodes.
Definition: pzgmesh.cpp:197
void Print(const TPZCompMesh &mesh, std::ostream &out=std::cout)
Print the information for the connect element.
Definition: pzconnect.cpp:67
REAL MaximumRadiusOfMesh()
Definition: pzcmesh.cpp:1833
void ConnectSolution(std::ostream &out)
Print the solution by connect index.
Definition: pzcmesh.cpp:1373
Defines the interface of a computational element. Computational Element.
Definition: pzcompel.h:59
bool IsCondensed() const
Access method to return the indication whether the connect is condensed or not.
Definition: pzconnect.h:223
virtual REAL CompareElement(int var, char *matname)
This method computes the norm of the difference of a post processed variable with @ the same post pro...
Definition: pzcompel.cpp:593
Contains the declaration of the TPZCondensedCompEl class, which implements an computational element w...
virtual int NSubElements() const =0
Returns the number of subelements of the element independent of the fact whether the element has alr...
TPZGeoEl * Father() const
Computes the normal vectors needed for forming HDiv vector valued shape functions.
Definition: pzgeoel.cpp:2574
int GetDefaultOrder()
Definition: pzcmesh.h:398
void TransferMultiphysicsSolution()
Transfer multiphysics mesh solution.
Definition: pzcmesh.cpp:470
Contains TPZShapeQuad class which implements the shape functions of a quadrilateral element...
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
Contains the TPZRefCube class which implements the uniform refinement of a geometric hexahedral eleme...
void Discontinuous2Continuous(int64_t disc_index, int64_t &new_index)
This method convert a discontinuous element with index disc_index in continuous element.
Definition: pzcmesh.cpp:1234
Contains the TPZRefPoint class which implements the uniform refinement of a geometric point element...
TPZFMatrix< STATE > & Solution()
Access the solution vector.
Definition: pzcmesh.h:219
static void WritePointer(const TPZSavable *obj, TPZStream *stream)
virtual void AutoBuild()
Creates the computational elements, and the degree of freedom nodes.
Definition: pzcmesh.h:474
Contains TPZShapeTriang class which implements the shape functions of a triangular element...
Contains TPZMetis class which implements the renumbering for elements of a mesh to minimize the band...
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzcmesh.cpp:1975
Implements computational element based on an interpolation space. Computational Element.
Definition: pzintel.h:27
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzblock.cpp:496
virtual MElementType Type() override
Type of the element.
TPZAdmChunkVector< TPZConnect > fConnectVec
List of pointers to nodes.
Definition: pzcmesh.h:62
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZMaterial * Material() const
Definition: pzbndcond.h:263
int Resequence(const int start=0)
Resequences blocks positioning.
Definition: pzblock.cpp:150
Implements an agglomerated discontinuous element. Computational Element.
virtual void Print(std::ostream &out=std::cout) const override
Prints the relevant data of the element to the output stream.
virtual void Print(std::ostream &out=std::cout) const
Prints mesh data.
Definition: pzcmesh.cpp:236
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
virtual int64_t GetFromSuperMesh(int64_t superind, TPZCompMesh *super)
Get an external connection from the supermesh - Supermesh is one mesh who contains the analised subme...
Definition: pzcmesh.cpp:1517
TPZAdmChunkVector< TPZGeoEl * > & ElementVec()
Methods for handling pzlists.
Definition: pzgmesh.h:138
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzadmchunk.h:132
void Permute(TPZVec< int64_t > &permute)
Permute the sequence number of the connect objects It is a permute gather operation.
Definition: pzcmesh.cpp:1321