NeoPZ
pzinterpolationspace.cpp
Go to the documentation of this file.
1 
6 #include "pzinterpolationspace.h"
7 #include "pzmaterialdata.h"
8 #include "pzbndcond.h"
9 #include "pzelmat.h"
10 #include "pzquad.h"
11 #include "TPZCompElDisc.h"
12 #include "TPZInterfaceEl.h"
13 #include "pztransfer.h"
14 #include "tpzchangeel.h"
15 
16 #include "pzlog.h"
17 
18 #ifdef LOG4CXX
19 static LoggerPtr logger(Logger::getLogger("pz.mesh.TPZInterpolationSpace"));
20 #endif
21 
22 TPZInterpolationSpace::TPZInterpolationSpace()
23 : TPZCompEl()
24 {
25  fPreferredOrder = -1;
26 }
27 
28 TPZInterpolationSpace::TPZInterpolationSpace(TPZCompMesh &mesh, const TPZInterpolationSpace &copy)
29 : TPZCompEl(mesh, copy)
30 {
31  fPreferredOrder = copy.fPreferredOrder;
32 }
33 
34 TPZInterpolationSpace::TPZInterpolationSpace(TPZCompMesh &mesh, const TPZInterpolationSpace &copy, std::map<int64_t,int64_t> &gl2lcElMap)
35 : TPZCompEl(mesh, copy, gl2lcElMap)
36 {
37  fPreferredOrder = copy.fPreferredOrder;
38 }
39 
40 TPZInterpolationSpace::TPZInterpolationSpace(TPZCompMesh &mesh, const TPZInterpolationSpace &copy, int64_t &index)
41 : TPZCompEl(mesh, copy, index)
42 {
43  fPreferredOrder = copy.fPreferredOrder;
44 }
45 
46 TPZInterpolationSpace::TPZInterpolationSpace(TPZCompMesh &mesh, TPZGeoEl *gel, int64_t &index)
47 : TPZCompEl(mesh,gel,index)
48 {
49  fPreferredOrder = mesh.GetDefaultOrder();
50 }
51 
52 TPZInterpolationSpace::~TPZInterpolationSpace(){
53 }
54 
55 int TPZInterpolationSpace::MaxOrder(){
56  const int n = this->NConnects();
57  int result = -1;
58  int side;
59  for(int i = 0; i < n; i++){
60  // skip unidentified connects
61  if (ConnectIndex(i) < 0) {
62  continue;
63  }
64  side = this->Connect(i).Order();
65  if (side > result) result = side;
66  }//i
67  return result;
68 }
69 
71 void TPZInterpolationSpace::AdjustIntegrationRule()
72 {
73  int order = MaxOrder();
74  int integrationruleorder = 0;
75  TPZMaterial * mat = this->Material();
76  if (mat) {
77  integrationruleorder = mat->IntegrationRuleOrder(order);
78  }else
79  {
80  integrationruleorder = order + order;
81  }
82  SetIntegrationRule(integrationruleorder);
83 }
84 
85 int TPZInterpolationSpace::ComputeIntegrationOrder() const {
86  DebugStop();
87  return 0;
88 }
89 
90 void TPZInterpolationSpace::Print(std::ostream &out) const {
91  out << __PRETTY_FUNCTION__ << std::endl;
92  TPZCompEl::Print(out);
93  out << "PreferredSideOrder " << fPreferredOrder << std::endl;
94 }
95 
96 void TPZInterpolationSpace::ComputeShape(TPZVec<REAL> &intpoint, TPZVec<REAL> &X,
97  TPZFMatrix<REAL> &jacobian, TPZFMatrix<REAL> &axes,
98  REAL &detjac, TPZFMatrix<REAL> &jacinv,
99  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi, TPZFMatrix<REAL> &dphidx){
100  TPZGeoEl * ref = this->Reference();
101  if (!ref){
102  PZError << "\nERROR AT " << __PRETTY_FUNCTION__ << " - this->Reference() == NULL\n";
103  return;
104  }//if
105 
106  ref->Jacobian( intpoint, jacobian, axes, detjac , jacinv);
107  this->Shape(intpoint,phi,dphi);
108  this->Convert2Axes(dphi, jacinv, dphidx);
109 
110 }
111 
112 void TPZInterpolationSpace::ComputeShape(TPZVec<REAL> &intpoint, TPZMaterialData &data){
113 
114 
115  this->ComputeShape(intpoint,data.x,data.jacobian,data.axes,data.detjac,data.jacinv,data.phi,data.dphi,data.dphix);
116 
117 }
118 
119 REAL TPZInterpolationSpace::InnerRadius(){
120  if (!this->Reference()){
121  PZError << "Error at " << __PRETTY_FUNCTION__ << " this->Reference() == NULL\n";
122  return 0.;
123  }
124  return this->Reference()->ElementRadius();
125 }
126 
127 void TPZInterpolationSpace::InitMaterialData(TPZMaterialData &data){
128  data.gelElId = this->Reference()->Id();
129  TPZMaterial *mat = Material();
130 #ifdef PZDEBUG
131  if(!mat)
132  {
133  mat= Material();
134  DebugStop();
135  }
136 #endif
137  mat->FillDataRequirements(data);
138  const int dim = this->Dimension();
139  const int nshape = this->NShapeF();
140  const int nstate = this->Material()->NStateVariables();
141  data.fShapeType = TPZMaterialData::EScalarShape;
142  data.phi.Redim(nshape,1);
143  data.dphi.Redim(dim,nshape);
144  data.dphix.Redim(dim,nshape);
145  data.axes.Redim(dim,3);
146  data.jacobian.Redim(dim,dim);
147  data.jacinv.Redim(dim,dim);
148  data.x.Resize(3);
149  if (data.fNeedsSol){
150  int64_t nsol = data.sol.size();
151  for (int64_t is=0; is<nsol; is++) {
152  data.sol[is].Resize(nstate);
153  data.dsol[is].Redim(dim,nstate);
154  }
155  }
156 }//void
157 
158 void TPZInterpolationSpace::ComputeRequiredData(TPZMaterialData &data,
159  TPZVec<REAL> &qsi){
160  data.intGlobPtIndex = -1;
161  this->ComputeShape(qsi, data);
162 
163  if (data.fNeedsSol){
164  if (data.phi.Rows()){//if shape functions are available
165  this->ComputeSolution(qsi, data);
166  }
167  else{//if shape functions are not available
168  this->ComputeSolution(qsi, data.sol, data.dsol, data.axes);
169  }
170  }//fNeedsSol
171 
172  data.x.Resize(3, 0.0);
173  Reference()->X(qsi, data.x);
174  data.xParametric = qsi;
175  TPZManVector<REAL,3> x_center(3,0.0);
176  int dim = Reference()->Dimension();
177  TPZVec<REAL> center_qsi(dim,0.0);
178 
179  if (Reference()->Dimension() == 2){
180  if (dim == EQuadrilateral) {
181  center_qsi[0] = 0.0;
182  center_qsi[1] = 0.0;
183  } else if (Reference()->Type() == ETriangle) {
184  center_qsi[0] = 0.25;
185  center_qsi[1] = 0.25;
186  }
187  }
188 
189  Reference()->X(center_qsi, x_center);
190  data.XCenter = x_center;
191 
192  if (data.fNeedsHSize){
193  data.HSize = 2.*this->InnerRadius();
194  }//fNeedHSize
195 
196  if (data.fNeedsNormal){
197  this->ComputeNormal(data);
198  }//fNeedsNormal
199 
200 }//void
201 
202 
203 void TPZInterpolationSpace::ComputeNormal(TPZMaterialData & data)
204 {
205  data.normal.Resize(3,0.);
206 
207  int thisFace, neighbourFace, i, dim;
208  TPZGeoEl * thisGeoEl, * neighbourGeoEl;
209 
210  thisGeoEl = this->Reference();
211  int thiseldim = thisGeoEl->Dimension();
212  TPZManVector<REAL,3> thisCenter(thiseldim,0.), thisXVol(3,0.), neighbourXVol(3,0.), vec(3), axes1(3), axes2(3);
213 
214  thisFace = thisGeoEl->NSides() - 1;
215  TPZGeoElSide thisside(thisGeoEl,thisFace);
216  TPZCompMesh *cmesh = this->Mesh();
217 
218 
219  TPZGeoElSide neighbourGeoElSide = thisGeoEl->Neighbour(thisFace);
220  int matid = neighbourGeoElSide.Element()->MaterialId();
221  while (neighbourGeoElSide != thisside) {
222  matid = neighbourGeoElSide.Element()->MaterialId();
223  if(cmesh->FindMaterial(matid) && neighbourGeoElSide.Element()->Dimension() > thiseldim)
224  {
225  break;
226  }
227  neighbourGeoElSide = neighbourGeoElSide.Neighbour();
228  }
229  neighbourGeoEl = neighbourGeoElSide.Element();
230  neighbourFace = neighbourGeoEl->NSides() - 1;
231 
232 
233  if(neighbourGeoEl == thisGeoEl)
234  {
235  // normal evaluation makes no sense since the internal element side doesn't have a neighbour.
236  return; // place a breakpoint here if this is an issue
237  }
238  int neightdim = neighbourGeoEl->Dimension();
239  TPZManVector<REAL,3> neighbourCenter(neightdim,0.);
240 
241  thisGeoEl->CenterPoint(thisFace, thisCenter);
242  neighbourGeoEl->CenterPoint(neighbourFace, neighbourCenter);
243 
244  thisGeoEl->X(thisCenter,thisXVol);
245  neighbourGeoEl->X(neighbourCenter,neighbourXVol);
246 
247  for(i = 0; i < 3; i++)
248  vec[i] = -neighbourXVol[i] + thisXVol[i];// vector towards the center of the neighbour element
249 
250  dim = thisGeoEl->Dimension();
251 
252  switch(dim)
253  {
254  case(0): // normal points towards the x-direction
255  data.normal[0] = 1.;
256  data.normal[1] = 0.;
257  data.normal[2] = 0.;
258  break;
259  case(1):
260  for(i = 0 ; i < 3; i ++) axes1[i] = data.axes(0,i); // rib direction
261  this->VectorialProd(axes1, vec, axes2);
262  this->VectorialProd(axes2, axes1, data.normal, true);
263  break;
264  case(2):
265  for(i = 0; i < 3; i++)
266  {
267  axes1[i] = data.axes(0,i);
268  axes2[i] = data.axes(1,i);
269  }
270  this->VectorialProd(axes1, axes2, data.normal, true);
271  break;
272  case(3):// in this case the normal becomes senseless. A null vector is passed instead
273  break;
274  default:
275  PZError << "TPZInterpolationSpace::ComputeNormal - unhandled element dimension\n";
276  }
277 
278  // ensuring the normal vector points towards the neighbour element
279 
280  REAL dot = 0.;
281  for(i = 0; i < 3; i++) dot += data.normal[i] * vec[i];
282 
283  if(dot < 0.)
284  for(i = 0; i < 3; i++) data.normal[i] *= -1.;
285 
286 }
287 
288 void TPZInterpolationSpace::VectorialProd(TPZVec<REAL> & ivec, TPZVec<REAL> & jvec, TPZVec<REAL> & kvec, bool unitary)
289 {
290  kvec.Resize(3);
291  kvec[0] = ivec[1]*jvec[2] - ivec[2]*jvec[1];
292  kvec[1] = -ivec[0]*jvec[2] + ivec[2]*jvec[0];
293  kvec[2] = ivec[0]*jvec[1] - ivec[1]*jvec[0];
294 
295  if(unitary)
296  {
297  REAL size = 0.;
298  int i;
299  for(i = 0; i < 3; i++)size += kvec[i] * kvec[i];
300  size = sqrt(size);
301  //if(size <= 1.e-9)PZError << "\nTPZInterpolationSpace::VectorialProd - null result\n";
302  for(i = 0; i < 3; i++)kvec[i] /= size;
303  }
304 }
305 
306 void TPZInterpolationSpace::CalcStiff(TPZElementMatrix &ek, TPZElementMatrix &ef){
307 
308  TPZMaterial * material = Material();
309  if(!material){
310  PZError << "Error at " << __PRETTY_FUNCTION__ << " this->Material() == NULL\n";
311  int matid = Reference()->MaterialId();
312  PZError << "Material Id which is missing = " << matid << std::endl;
313  ek.Reset();
314  ef.Reset();
315  return;
316  }
317 
318 #ifdef LOG4CXX
319  if (logger->isDebugEnabled())
320  {
321  std::stringstream sout;
322  sout << __PRETTY_FUNCTION__ << " material id " << material->Id();
323  LOGPZ_DEBUG(logger,sout.str());
324  }
325 #endif
326  this->InitializeElementMatrix(ek,ef);
327 
328  if (this->NConnects() == 0) return;//boundary discontinuous elements have this characteristic
329 
330  TPZMaterialData data;
331  this->InitMaterialData(data);
332  data.p = this->MaxOrder();
333 
334  int dim = Dimension();
335  TPZManVector<REAL,3> intpoint(dim,0.);
336  REAL weight = 0.;
337 
338  TPZAutoPointer<TPZIntPoints> intrule = GetIntegrationRule().Clone();
339 
340 //#ifdef PZDEBUG
341 // {
342 
343 // TPZManVector<int,3> intorder(dim,this->MaxOrder()*2);
344 // TPZAutoPointer<TPZIntPoints> intrule_clone = intrule->Clone();
345 // intrule_clone->SetOrder(intorder);
346 //
347 // if(intrule_clone->NPoints() > intrule->NPoints()){
348 // std::cout << "Element " << fIndex << " Bad integration rule points needed = " << intrule_clone->NPoints() << "; points obtained = " << intrule->NPoints() << std::endl;
349 // intrule = intrule_clone;
350 // }
351 
352 // }
353 //#endif
354 
355 // if(material->HasForcingFunction())
356 // {
357 // int maxorder = intrule->GetMaxOrder();
358 // if (maxorder > order) {
359 // order = maxorder;
360 // }
361 // }
362 // TPZManVector<int,3> intorder(dim,order);
363 // intrule->SetOrder(intorder);
364 
365  int intrulepoints = intrule->NPoints();
366  for(int int_ind = 0; int_ind < intrulepoints; ++int_ind){
367  intrule->Point(int_ind,intpoint,weight);
368  data.intLocPtIndex = int_ind;
369  this->ComputeRequiredData(data, intpoint);
370  weight *= fabs(data.detjac);
371 
372  material->Contribute(data, weight, ek.fMat, ef.fMat);
373  }//loop over integratin points
374 
375 }//CalcStiff
376 
377 void TPZInterpolationSpace::CalcResidual(TPZElementMatrix &ef){
378 
379  TPZMaterial * material = Material();
380  if(!material){
381  PZError << "Error at " << __PRETTY_FUNCTION__ << " this->Material() == NULL\n";
382  ef.Reset();
383  return;
384  }
385 
386  this->InitializeElementMatrix(ef);
387 
388  if (this->NConnects() == 0) return; //boundary discontinuous elements have this characteristic
389 
390  TPZMaterialData data;
391  this->InitMaterialData(data);
392  data.p = this->MaxOrder();
393 
394  int dim = Dimension();
395  TPZManVector<REAL,3> intpoint(dim,0.);
396  REAL weight = 0.;
397 
398  TPZAutoPointer<TPZIntPoints> intrule = GetIntegrationRule().Clone();
399 // if(material->HasForcingFunction())
400 // {
401 // int maxorder = intrule->GetMaxOrder();
402 // if (maxorder > order) {
403 // order = maxorder;
404 // }
405 // }
406 // TPZManVector<int,3> intorder(dim,order);
407 // intrule->SetOrder(intorder);
408  // material->SetIntegrationRule(intrule, data.p, dim);
409 
410  int intrulepoints = intrule->NPoints();
411  for(int int_ind = 0; int_ind < intrulepoints; ++int_ind){
412  intrule->Point(int_ind,intpoint,weight);
413 
414  //this->ComputeShape(intpoint, data.x, data.jacobian, data.axes, data.detjac, data.jacinv, data.phi, data.dphix);
415 
416  //weight *= fabs(data.detjac);
417 
418  data.intLocPtIndex = int_ind;
419 
420  this->ComputeRequiredData(data, intpoint);
421  weight *= fabs(data.detjac);
422 
423  material->Contribute(data,weight,ef.fMat);
424 
425  }//loop over integratin points
426 
427 }//CalcResidual
428 
429 void TPZInterpolationSpace::InitializeElementMatrix(TPZElementMatrix &ek, TPZElementMatrix &ef){
430  TPZMaterial *mat = this->Material();
431  const int numdof = mat->NStateVariables();
432  const int nshape = this->NShapeF();
433  const int ncon = this->NConnects();
434  const int numloadcases = mat->NumLoadCases();
435 
436  ek.fMesh = Mesh();
437  ek.fType = TPZElementMatrix::EK;
438  ek.fOneRestraints = GetShapeRestraints();
439  ef.fOneRestraints = ek.fOneRestraints;
440 
441  ef.fMesh = Mesh();
442  ef.fType = TPZElementMatrix::EF;
443 
444  ek.fBlock.SetNBlocks(ncon);
445  ef.fBlock.SetNBlocks(ncon);
446 
447  int i;
448  int numeq=0;
449  for(i=0; i<ncon; i++){
450  TPZConnect &c = Connect(i);
451  int nshape = c.NShape();
452 #ifdef PZDEBUG
453  if (nshape != NConnectShapeF(i,c.Order())) {
454  DebugStop();
455  }
456 #endif
457  int nstate = c.NState();
458 
459 #ifdef PZDEBUG
460  if(nstate != numdof)
461  {
462  DebugStop();
463  }
464 #endif
465  ek.fBlock.Set(i,nshape*nstate);
466  ef.fBlock.Set(i,nshape*nstate);
467  numeq += nshape*nstate;
468  }
469  ek.fMat.Redim(numeq,numeq);
470  ef.fMat.Redim(numeq,numloadcases);
471  ek.fConnect.Resize(ncon);
472  ef.fConnect.Resize(ncon);
473  for(i=0; i<ncon; i++){
474  (ef.fConnect)[i] = ConnectIndex(i);
475  (ek.fConnect)[i] = ConnectIndex(i);
476  }
477 }//void
478 
479 void TPZInterpolationSpace::InitializeElementMatrix(TPZElementMatrix &ef){
480  TPZMaterial *mat = this->Material();
481  const int numdof = mat->NStateVariables();
482  const int ncon = this->NConnects();
483  const int nshape = this->NShapeF();
484  const int numeq = nshape*numdof;
485  const int numloadcases = mat->NumLoadCases();
486  ef.fMesh = Mesh();
487  ef.fType = TPZElementMatrix::EF;
488  ef.fMat.Redim(numeq,numloadcases);
489  ef.fBlock.SetNBlocks(ncon);
490  ef.fOneRestraints = GetShapeRestraints();
491  int i;
492  for(i=0; i<ncon; i++){
493  TPZConnect &c = Connect(i);
494  unsigned int nshapec = c.NShape();
495 #ifdef PZDEBUG
496  if (NConnectShapeF(i, c.Order()) != nshapec || c.NState() != numdof) {
497  DebugStop();
498  }
499 #endif
500  ef.fBlock.Set(i,nshapec*numdof);
501  }
502  ef.fConnect.Resize(ncon);
503  for(i=0; i<ncon; i++){
504  (ef.fConnect)[i] = ConnectIndex(i);
505  }
506 }//void
507 
508 void TPZInterpolationSpace::Solution(TPZVec<REAL> &qsi,int var,TPZVec<STATE> &sol) {
509  if(var >= 100) {
510  TPZCompEl::Solution(qsi,var,sol);
511  return;
512  }
513  if(var == 99) {
514  sol[0] = GetPreferredOrder();
515  // if (sol[0] != 2) {
516  // std::cout << __PRETTY_FUNCTION__ << " preferred order " << sol[0] << std::endl;
517  // }
518  return;
519  }
520 
521  TPZMaterial * material = this->Material();
522  if(!material) {
523  sol.Resize(0);
524  return;
525  }
526 
527  TPZMaterialData data;
528  this->InitMaterialData(data);
529  data.p = this->MaxOrder();
530  this->ComputeShape(qsi, data);
531  this->ComputeSolution(qsi,data);
532 
533  data.x.Resize(3);
534  this->Reference()->X(qsi, data.x);
535 
536  int solSize = material->NSolutionVariables(var);
537  sol.Resize(solSize);
538  sol.Fill(0.);
539  material->Solution(data, var, sol);
540 }
541 
542 void TPZInterpolationSpace::InterpolateSolution(TPZInterpolationSpace &coarsel){
543  // accumulates the transfer coefficients between the current element and the
544  // coarse element into the transfer matrix, using the transformation t
545  TPZMaterial * material = Material();
546  if(!material) {
547  PZError << __PRETTY_FUNCTION__ << " this->Material() == NULL " << std::endl;
548  return;
549  }
550 
551  TPZTransform<> t(Dimension());
552  TPZGeoEl *ref = Reference();
553 
554  //Cedric 16/03/99
555  // Reference()->BuildTransform(NConnects(),coarsel.Reference(),t);
556  t = Reference()->BuildTransform2(ref->NSides()-1,coarsel.Reference(),t);
557 
558  int locnod = NConnects();
559  // int cornod = coarsel.NConnects();
560  int locmatsize = NShapeF();
561  int cormatsize = coarsel.NShapeF();
562  int nvar = material->NStateVariables();
563  int dimension = Dimension();
564  if (!dimension) {
565  PZError << "\nExiting " << __PRETTY_FUNCTION__ << " - trying to interpolate a node solution.\n";
566  return ;
567  }
568 
569  //TPZFMatrix<REAL> loclocmat(locmatsize,locmatsize,0.);
570  TPZFMatrix<STATE> loclocmat(locmatsize,locmatsize,0.);
571  //TPZFMatrix<REAL> projectmat(locmatsize,nvar,0.);
572  TPZFMatrix<STATE> projectmat(locmatsize,nvar,0.);
573 
574  TPZManVector<int,3> prevorder(dimension);
575  TPZAutoPointer<TPZIntPoints> intrule = GetIntegrationRule().Clone();
576  intrule->GetOrder(prevorder);
577 
578  int thismaxorder = this->MaxOrder();
579  int coarsemaxorder = coarsel.MaxOrder();
580  int maxorder = (thismaxorder > coarsemaxorder) ? thismaxorder : coarsemaxorder;
581  // Cesar 2003-11-25 -->> To avoid integration warnings...
582  maxorder = (2*maxorder > intrule->GetMaxOrder() ) ? intrule->GetMaxOrder() : 2*maxorder;
583  TPZManVector<int,3> order(dimension,maxorder);
584  for(int dim = 0; dim < dimension; dim++) {
585  order[dim] = maxorder*2;
586  }
587  intrule->SetOrder(order);
588 
589  TPZFNMatrix<220> locphi(locmatsize,1);
590  TPZFNMatrix<660> locdphi(dimension,locmatsize);//derivative of the shape function in the master domain
591  TPZFNMatrix<660> locdphidx(dimension,locmatsize);//derivative of the shape function in the deformed domain
592 
593  TPZFNMatrix<220> corphi(cormatsize,1);
594  TPZFNMatrix<660> cordphi(dimension,cormatsize);//derivative of the shape function in the master domain
595 
596  TPZManVector<REAL,3> int_point(dimension),coarse_int_point(dimension);
597  TPZFNMatrix<9> jacobian(dimension,dimension),jacinv(dimension,dimension);
598  TPZFNMatrix<9> axes(3,3,0.), coarseaxes(3,3,0.);
599  REAL zero = 0.;
600  TPZManVector<REAL,3> x(3,zero);
601  //TPZManVector<TPZManVector<REAL,10>, 10> u(1);
602  TPZSolVec u(1);
603  u[0].resize(nvar);
604  //TPZManVector<TPZFNMatrix<30>, 10> du(1);
605  TPZGradSolVec du(1);
606  du[0].Redim(dimension,nvar);
607 
608  int numintpoints = intrule->NPoints();
609  REAL weight;
610  int lin,ljn,cjn;
611  TPZConnect *df;
612  // TPZBlock &coarseblock = coarsel.Mesh()->Block();
613 
614  for(int int_ind = 0; int_ind < numintpoints; ++int_ind) {
615  intrule->Point(int_ind,int_point,weight);
616  REAL jac_det = 1.;
617  this->ComputeShape(int_point, x, jacobian, axes, jac_det, jacinv, locphi, locdphi, locdphidx);
618  weight *= jac_det;
619  t.Apply(int_point,coarse_int_point);
620  coarsel.ComputeSolution(coarse_int_point, u, du, coarseaxes);
621 
622  for(lin=0; lin<locmatsize; lin++) {
623  for(ljn=0; ljn<locmatsize; ljn++) {
624  loclocmat(lin,ljn) += weight*locphi(lin,0)*locphi(ljn,0);
625  }
626  for(cjn=0; cjn<nvar; cjn++) {
627  projectmat(lin,cjn) += (STATE)weight*(STATE)locphi(lin,0)*u[0][cjn];
628  }
629  }
630  jacobian.Zero();
631  }
632  REAL large = 0.;
633  for (lin = 0; lin < locmatsize; lin++) {
634  for (cjn = 0; cjn < locmatsize; cjn++) {
635  REAL val = fabs(loclocmat(lin,cjn));
636  large = large<val ? val : large;
637  }
638  }
639  loclocmat *= 1./large;
640  projectmat *= 1./large;
641 
642  loclocmat.SolveDirect(projectmat,ELU);
643  // identify the non-zero blocks for each row
644  TPZBlock<STATE> &fineblock = Mesh()->Block();
645  int iv=0,in;
646  for(in=0; in<locnod; in++) {
647  df = &Connect(in);
648  int64_t dfseq = df->SequenceNumber();
649  int dfvar = fineblock.Size(dfseq);
650  for(ljn=0; ljn<dfvar; ljn++) {
651  fineblock(dfseq,0,ljn,0) = projectmat(iv/nvar,iv%nvar);
652  iv++;
653  }
654  }
655  intrule->SetOrder(prevorder);
656 }//InterpolateSolution
657 
658 void TPZInterpolationSpace::CreateInterfaces(bool BetweenContinuous) {
659  //nao verifica-se caso o elemento de contorno
660  //eh maior em tamanho que o interface associado
661  //caso AdjustBoundaryElement nao for aplicado
662  //a malha eh criada consistentemente
663  TPZGeoEl *ref = Reference();
664  int nsides = ref->NSides();
665  int InterfaceDimension = this->Material()->Dimension() - 1;
666  int side;
667  nsides--;//last face
668  for(side=nsides;side>=0;side--) {
669  if(ref->SideDimension(side) != InterfaceDimension) continue;
670  TPZCompElSide thisside(this,side);
671  if(this->ExistsInterface(thisside.Reference())) {
672  // cout << "TPZCompElDisc::CreateInterface inconsistent: interface already exists\n";
673  continue;
674  }
675  TPZStack<TPZCompElSide> highlist;
676  thisside.HigherLevelElementList(highlist,0,1);
677  //a interface se cria uma vez so quando existem ambos
678  //elementos esquerdo e direito (compu tacionais)
679  if(!highlist.NElements()) {
680  this->CreateInterface(side, BetweenContinuous);//s�tem iguais ou grande => pode criar a interface
681  } else {
682  int64_t ns = highlist.NElements();
683  int64_t is;
684  for(is=0; is<ns; is++) {//existem pequenos ligados ao lado atual
685  const int higheldim = highlist[is].Reference().Dimension();
686  if(higheldim != InterfaceDimension) continue;
687  // TPZCompElDisc *del = dynamic_cast<TPZCompElDisc *> (highlist[is].Element());
688  // if(!del) continue;
689 
690  TPZCompEl *del = highlist[is].Element();
691  if(!del) continue;
692 
693  TPZCompElSide delside( del, highlist[is].Side() );
694  TPZInterpolationSpace * delSp = dynamic_cast<TPZInterpolationSpace*>(del);
695  if (!delSp) {
696  PZError << "\nERROR AT " << __PRETTY_FUNCTION__ << " - CASE NOT AVAILABLE\n";
697  return;
698  }
699  if ( delSp->ExistsInterface(delside.Reference()) ) {
700  // cout << "TPZCompElDisc::CreateInterface inconsistent: interface already exists\n";
701  }
702  else {
703  delSp->CreateInterface(highlist[is].Side(), BetweenContinuous);
704  }
705  }
706  }
707  }
708 }
709 
710 TPZInterfaceElement * TPZInterpolationSpace::CreateInterface(int side, bool BetweenContinuous)
711 {
712  // LOGPZ_INFO(logger, "Entering CreateInterface");
713 
714  TPZGeoEl *ref = Reference();
715  if(!ref) {
716  LOGPZ_WARN(logger, "Exiting CreateInterface Null reference reached - NULL interface returned");
717  return NULL;
718  }
719 
720  TPZCompElSide thisside(this,side);
721  TPZStack<TPZCompElSide> list;
722  list.Resize(0);
723  thisside.EqualLevelElementList(list,0,1);//retorna distinto ao atual ou nulo
724  const int64_t size = list.NElements();
725 
726  if (size > 1) {
727  DebugStop();
728  }
729  //espera-se ter os elementos computacionais esquerdo e direito
730  //ja criados antes de criar o elemento interface
731  if(size){
732  //Interface has the same material of the neighbour with lesser dimension.
733  //It makes the interface have the same material of boundary conditions (TPZCompElDisc with interface dimension)
734  int64_t index;
735 
736  TPZCompEl *list0 = list[0].Element();
737  int list0side = list[0].Side();
738  TPZCompElDisc * thisdisc = dynamic_cast<TPZCompElDisc*>(this);
739  TPZCompElDisc * neighdisc = dynamic_cast<TPZCompElDisc*>(list0);
740  int thisside = side;
741  int neighside = list0side;
742 
743  if (BetweenContinuous == false){
744  //It means at least one element must be discontinuous
745  if (!thisdisc && !neighdisc){
746  return NULL;
747  }
748  }
749  TPZInterfaceElement * newcreatedinterface = NULL;
750 
751  if (Dimension() == list0->Dimension())
752  {
753  const int matid = this->Mesh()->Reference()->InterfaceMaterial(this->Material()->Id(), list0->Material()->Id() );
754  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid); //isto acertou as vizinhanas da interface geometrica com o atual
755  if(!gel)
756  {
757 #ifdef LOG4CXX
758  if (logger->isDebugEnabled())
759  {
760  std::stringstream sout;
761  sout << "CreateBCGeoEl devolveu zero!@@@@";
762  LOGPZ_DEBUG(logger,sout.str());
763  }
764 #endif
765  DebugStop();
766  }
767  TPZCompElSide thiscompelside(this, thisside);
768  TPZCompElSide neighcompelside(list0, neighside);
769  int64_t index;
770  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,thiscompelside,neighcompelside);
771 
772  } else
773  {
774 
775  if(Dimension() > list0->Dimension())
776  {
777  const int matid = list0->Material()->Id();
778  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid); //isto acertou as vizinhanas da interface geometrica com o atual
779  if (!gel) {
780  DebugStop();
781  }
782 
783  //o de volume eh o direito caso um deles seja BC
784  //a normal aponta para fora do contorno
785  TPZCompElSide thiscompelside(this, thisside);
786  TPZCompElSide neighcompelside(list0, neighside);
787  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,thiscompelside,neighcompelside);
788  } else {
789  const int matid = this->Material()->Id();
790  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid); //isto acertou as vizinhanas da interface geometrica com o atual
791  if(!gel) DebugStop();
792  //caso contrario ou caso ambos sejam de volume
793  TPZCompElSide thiscompelside(this, thisside);
794  TPZCompElSide neighcompelside(list0, neighside);
795  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,neighcompelside,thiscompelside);
796  }
797  }
798 
799 
801 #ifdef PZDEBUG
802  {
803  TPZGeoEl * faceGel = newcreatedinterface->Reference();
804  TPZGeoEl * leftGel = newcreatedinterface->LeftElement()->Reference();
805  TPZGeoEl * rightGel = newcreatedinterface->RightElement()->Reference();
806  bool leftIsLinear = leftGel->IsLinearMapping();
807  bool rightIsLinear = rightGel->IsLinearMapping();
808  if(!leftIsLinear && !rightIsLinear){
809  if(faceGel->IsGeoBlendEl() == false){
810  std::cout << "\nError at " << __PRETTY_FUNCTION__ << "\n";
811 #ifdef LOG4CXX
812  {
813  std::stringstream sout;
814  sout << "\nError at " << __PRETTY_FUNCTION__ << "\n";
815  sout << "left gel:\n";
816  leftGel->Print(sout);
817  sout << "right gel:";
818  rightGel->Print(sout);
819  sout << "face gel:";
820  faceGel->Print(sout);
821  LOGPZ_ERROR(logger,sout.str());
822  }
823 #endif
824  }
825  }
826  }
827 #endif
828 
830  return newcreatedinterface;
831  }
832 
833  //If there is no equal level element, we try the lower elements.
834  //Higher elements will not be considered by this method. In that case the interface must be created by the neighbour.
835  TPZCompElSide lower = thisside.LowerLevelElementList(0);
836  if(lower.Exists()){
837  //Interface has the same material of the neighbour with lesser dimension.
838  //It makes the interface has the same material of boundary conditions (TPZCompElDisc with interface dimension)
839  int matid;
840  int thisdim = this->Dimension();
841  int neighbourdim = lower.Element()->Dimension();
842 
843  if (thisdim == neighbourdim){
844  // matid = this->Material()->Id();
845  matid = this->Mesh()->Reference()->InterfaceMaterial(this->Material()->Id(), lower.Element()->Material()->Id() );
846  }
847  else { //one element is a boundary condition
848  if (thisdim < neighbourdim) matid = this->Material()->Id();
849  else
850  {
851  TPZCompEl *cel = lower.Element();
852  if (!cel) {
853  DebugStop();
854  }
855  TPZMaterial *mat = cel->Material();
856  if (!mat ) {
857  DebugStop();
858  }
859  matid = lower.Element()->Material()->Id();
860  }
861  }
862 
863 
864  TPZCompEl *lowcel = lower.Element();
865  int lowside = lower.Side();
866  TPZCompElDisc * thisdisc = dynamic_cast<TPZCompElDisc*>(this);
867  TPZCompElDisc * neighdisc = dynamic_cast<TPZCompElDisc*>(lowcel);
868  int thisside = side;
869  int neighside = lowside;
870 
871  if (BetweenContinuous == false){
872  //It means at least one element must be discontinuous
873  if (!thisdisc && !neighdisc){
874  return NULL;
875  }
876  }
877  TPZInterfaceElement * newcreatedinterface = NULL;
878 
879  int64_t index;
880 
881  if(Dimension() == lowcel->Dimension()){
882 
883  const int matid = this->Mesh()->Reference()->InterfaceMaterial(lowcel->Material()->Id(), this->Material()->Id() );
884  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid); //isto acertou as vizinhanas da interface geometrica com o atual
885  if(!gel) DebugStop();
886 
887  TPZCompElSide lowcelcompelside(lowcel, neighside);
888  TPZCompElSide thiscompelside(this, thisside);
889  int64_t index;
890  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,lowcelcompelside,thiscompelside);
891  }
892  else{
893 
894  if(Dimension() > lowcel->Dimension()){
895  const int matid = lowcel->Material()->Id();
896  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid);
897  if (!gel) {
898  DebugStop();
899  }
900 
901  //para que o elemento esquerdo seja de volume
902  TPZCompElSide thiscompelside(this, thisside);
903  TPZCompElSide lowcelcompelside(lowcel, neighside);
904  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,thiscompelside,lowcelcompelside);
905  } else {
906  const int matid = this->Material()->Id();
907  TPZGeoEl *gel = ref->CreateBCGeoEl(side,matid);
908  if (!gel) {
909  DebugStop();
910  }
911  TPZCompElSide thiscompelside(this, thisside);
912  TPZCompElSide lowcelcompelside(lowcel, neighside);
913 #ifdef LOG4CXX_KEEP
914  if (logger->isDebugEnabled())
915  {
916  std::stringstream sout;
917  sout << __PRETTY_FUNCTION__ << " left element";
918  sout << lowcelcompelside << thiscompelside;
919  sout << "Left Element ";
920  lowcelcompelside.Element()->Print(sout);
921  sout << "Right Element ";
922  thiscompelside.Element()->Print(sout);
923  LOGPZ_DEBUG(logger,sout.str())
924  }
925 #endif
926  newcreatedinterface = new TPZInterfaceElement(*fMesh,gel,index,lowcelcompelside,thiscompelside);
927  }
928  }
930 #ifdef PZDEBUG
931  {
932  TPZGeoEl * faceGel = newcreatedinterface->Reference();
933  TPZGeoEl * leftGel = newcreatedinterface->LeftElement()->Reference();
934  TPZGeoEl * rightGel = newcreatedinterface->RightElement()->Reference();
935  bool leftIsLinear = leftGel->IsLinearMapping();
936  bool rightIsLinear = rightGel->IsLinearMapping();
937  if(!leftIsLinear && !rightIsLinear){
938  if(faceGel->IsGeoBlendEl() == false){
939  std::cout << "\nError at " << __PRETTY_FUNCTION__ << "\n";
940 #ifdef LOG4CXX
941  {
942  std::stringstream sout;
943  sout << "\nError at " << __PRETTY_FUNCTION__ << "\n";
944  sout << "left gel:\n";
945  leftGel->Print(sout);
946  sout << "right gel:";
947  rightGel->Print(sout);
948  sout << "face gel:";
949  faceGel->Print(sout);
950  LOGPZ_ERROR(logger,sout.str());
951  }
952 #endif
953  }
954  }
955  }
956 #endif
957 
959  return newcreatedinterface;
960  }
961  return NULL;
962 }
963 
964 int TPZInterpolationSpace::ExistsInterface(TPZGeoElSide geosd){
965 
966  TPZGeoElSide neighside = geosd.Neighbour();
967  while(neighside.Element() && neighside.Element() != geosd.Element()){
968  TPZCompElSide neighcompside = neighside.Reference();
969  neighside = neighside.Neighbour();
970  if(!neighcompside.Element()) continue;
971  if(neighcompside.Element()->Type() == EInterface)
972  return 1;
973  }
974  return 0;
975 }
976 
977 void TPZInterpolationSpace::RemoveInterfaces(){
978 
979  int nsides = Reference()->NSides();
980  if (!this->Material()){
981  std::stringstream mess;
982  mess << __PRETTY_FUNCTION__ << " - this->Material() == NULL, I can't RemoveInterfaces()";
983  PZError << mess.str() << std::endl;
984  LOGPZ_ERROR(logger, mess.str());
985  return;
986  }
987  int InterfaceDimension = this->Material()->Dimension() - 1;
988  int is;
989  TPZStack<TPZCompElSide> list,equal;
990  for(is=0;is<nsides;is++){
991  TPZCompElSide thisside(this,is);
992  if(thisside.Reference().Dimension() != InterfaceDimension) continue;
993  // procurar na lista de elementos iguais
994  list.Resize(0);// o lado atual �uma face
995  //thisside.EqualLevelElementList(list,0,0);// monta a lista de elementos iguais
996  RemoveInterface(is);// chame remove interface do elemento atual (para o side atual)
997  thisside.HigherLevelElementList(list,0,0);// procurar na lista de elementos menores (todos)
998  int64_t size = list.NElements(), i; // 'isto pode incluir elementos interfaces'
999  //tirando os elementos de interface da lista
1000  for(i=0;i<size;i++){
1001  if(list[i].Element()->Type() == EInterface) {
1002  LOGPZ_DEBUG(logger, "Removing interface element from the list of higher level elements");
1003  //This need to be done because otherwise list could be invalidated when an interface is removed.
1004  list[i] = TPZCompElSide();//tirando interface
1005  }
1006  }
1007  for(i=0;i<size;i++){// percorre os elementos menores
1008  if(!list[i].Element()) continue;
1009  TPZGeoElSide geolist = list[i].Reference();//TESTE
1010  if(geolist.Dimension() != InterfaceDimension) continue;
1011  equal.Resize(0);// para cada elemento menor e' preciso verificar a dimensao,
1012  list[i].EqualLevelElementList(equal,0,0);//montar a lista de elementos iguais (todos)
1013  equal.Push(list[i]);//n� �incorporado no m�odo anterior
1014  int64_t neq = equal.NElements(),k=-1;
1015  while(++k < neq) if(equal[k].Element()->Type() != EInterface) break;//procurando elemento descont�uo cujo
1016  if(!neq || k == neq){ //lado faz parte da parti� do lado side do this
1017  LOGPZ_FATAL(logger, " Inconsistency of data");
1018  DebugStop();//elemento descont�uo n� achado: ERRO
1019  }// chame removeinterface do elemento menor
1020 
1021  TPZInterpolationSpace * equalkSp = dynamic_cast<TPZInterpolationSpace*>(equal[k].Element());
1022  if (!equalkSp){
1023  PZError << "\nERROR AT " << __PRETTY_FUNCTION__ << " - CASE NOT AVAILABLE\n";
1024  return;
1025  }
1026  equalkSp->RemoveInterface(equal[k].Side());
1027  }
1028  }
1029 
1030 }
1031 
1032 void TPZInterpolationSpace::RemoveInterface(int side) {
1033 
1034  TPZStack<TPZCompElSide> list;
1035  list.Resize(0);
1036  TPZCompElSide thisside(this,side);
1037  thisside.EqualLevelElementList(list,0,0);// monta a lista de elementos iguais
1038  int64_t size = list.NElements(),i=-1;
1039  while(++i < size) if(list[i].Element()->Type() == EInterface) break;// procura aquele que e derivado de TPZInterfaceEl
1040  if(!size || i == size){
1041 #ifdef LOG4CXX_keep
1042  if (logger->isDebugEnabled())
1043  {
1044  std::stringstream sout;
1045  sout << __PRETTY_FUNCTION__ << " no interface element found\n";
1046  Print(sout);
1047  LOGPZ_DEBUG(logger,sout.str())
1048  }
1049 #endif
1050  return;// nada a ser feito
1051  }
1052  // aqui existe a interface
1053  TPZCompEl *cel = list[i].Element();
1054 #ifdef LOG4CXX
1055  TPZGeoEl *gel = cel->Reference();
1056  if (logger->isDebugEnabled())
1057  {
1058  std::stringstream sout;
1059  sout << __PRETTY_FUNCTION__ << " element index " << Index() << " side " << std::endl;
1060  sout << "geometric element reference count " << gel->NumInterfaces();
1061  LOGPZ_DEBUG(logger,sout.str())
1062  }
1063 #endif
1064  delete cel;
1065 }
1066 
1067 void TPZInterpolationSpace::EvaluateError(std::function<void(const TPZVec<REAL> &loc,TPZVec<STATE> &val,TPZFMatrix<STATE> &deriv)> fp,
1068  TPZVec<REAL> &errors,bool store_error){
1069  TPZMaterial * material = Material();
1070  //TPZMaterial * matptr = material.operator->();
1071  if (!material) {
1072  PZError << "TPZInterpolatedElement::EvaluateError : no material for this element\n";
1073  Print(PZError);
1074  return;
1075  }
1076  if (dynamic_cast<TPZBndCond *>(material)) {
1077  LOGPZ_INFO(logger, "Exiting EvaluateError - null error - boundary condition material.");
1078  return;
1079  }
1080  int NErrors = material->NEvalErrors();
1081  errors.Resize(NErrors);
1082  errors.Fill(0.);
1083  int problemdimension = Mesh()->Dimension();
1084  if(Reference()->Dimension() < problemdimension) return;
1085 
1086  // Adjust the order of the integration rule
1087  int dim = Dimension();
1088  TPZAutoPointer<TPZIntPoints> intrule = this->GetIntegrationRule().Clone();
1089  int maxIntOrder = intrule->GetMaxOrder();
1090  TPZManVector<int,3> prevorder(dim), maxorder(dim, maxIntOrder);
1091  //end
1092  intrule->GetOrder(prevorder);
1093  const int order_limit = 8;
1094  if(maxIntOrder > order_limit)
1095  {
1096  if (prevorder[0] > order_limit) {
1097  maxIntOrder = prevorder[0];
1098  }
1099  else
1100  {
1101  maxIntOrder = order_limit;
1102  }
1103  }
1104 
1105  intrule->SetOrder(maxorder);
1106 
1107  int ndof = material->NStateVariables();
1108  int nflux = material->NFluxes();
1109  TPZManVector<STATE,10> u_exact(ndof);
1110  TPZFNMatrix<9,STATE> du_exact(dim+1,ndof);
1111  TPZManVector<REAL,10> intpoint(problemdimension), values(NErrors);
1112  REAL weight;
1113  TPZManVector<STATE,9> flux_el(nflux,0.);
1114 
1115  TPZMaterialData data;
1116  this->InitMaterialData(data);
1117  int nintpoints = intrule->NPoints();
1118 
1119  for(int nint = 0; nint < nintpoints; nint++) {
1120 
1121  values.Fill(0.0);
1122  intrule->Point(nint,intpoint,weight);
1123 
1124  //in the case of the hdiv functions
1125  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
1126  if(shapetype==data.EVecandShape){
1127  this->ComputeRequiredData(data, intpoint);
1128  }
1129  else
1130  {
1131  this->ComputeShape(intpoint, data.x, data.jacobian, data.axes, data.detjac, data.jacinv, data.phi, data.dphi, data.dphix);
1132  }
1133  weight *= fabs(data.detjac);
1134  // this->ComputeSolution(intpoint, data.phi, data.dphix, data.axes, data.sol, data.dsol);
1135  //this->ComputeSolution(intpoint, data);
1136  //contribuicoes dos erros
1137  TPZGeoEl * ref = this->Reference();
1138  ref->X(intpoint, data.x);
1139  if(fp) {
1140  fp(data.x,u_exact,du_exact);
1141 
1142  if(data.fVecShapeIndex.NElements())
1143  {
1144  this->ComputeSolution(intpoint, data);
1145 
1146  material->ErrorsHdiv(data,u_exact,du_exact,values);
1147 
1148  }
1149  else{
1150  this->ComputeSolution(intpoint, data.phi, data.dphix, data.axes, data.sol, data.dsol);
1151  material->Errors(data.x,data.sol[0],data.dsol[0],data.axes,flux_el,u_exact,du_exact,values);
1152  }
1153 
1154  for(int ier = 0; ier < NErrors; ier++)
1155  errors[ier] += weight*values[ier];
1156  }
1157 
1158  }//fim for : integration rule
1159 
1160  //Norma sobre o elemento
1161  for(int ier = 0; ier < NErrors; ier++){
1162  errors[ier] = sqrt(errors[ier]);
1163  }//for ier
1164  if(store_error)
1165  {
1166  int64_t index = Index();
1167  TPZFMatrix<STATE> &elvals = Mesh()->ElementSolution();
1168  if (elvals.Cols() < NErrors) {
1169  std::cout << "The element solution of the mesh should be resized before EvaluateError\n";
1170  DebugStop();
1171  }
1172  for (int ier=0; ier <NErrors; ier++) {
1173  elvals(index,ier) = errors[ier];
1174  }
1175  }
1176  intrule->SetOrder(prevorder);
1177 
1178 }//method
1179 
1180 void TPZInterpolationSpace::ComputeError(int errorid,
1181  TPZVec<REAL> &error){
1182 
1183  TPZMaterial * material = Material();
1184  if(!material){
1185  std::cout << "TPZCompElDisc::ComputeError : no material for this element\n";
1186  return;
1187  }
1188 
1189  TPZMaterialData data;
1190  this->InitMaterialData(data);
1191 
1192  REAL weight;
1193  int dim = Dimension();
1194  TPZVec<REAL> intpoint(dim,0.);
1195 
1196  TPZAutoPointer<TPZIntPoints> intrule = this->GetIntegrationRule().Clone();
1197 
1198  TPZManVector<int,3> prevorder(dim), maxorder(dim, this->MaxOrder());
1199  intrule->GetOrder(prevorder);
1200  intrule->SetOrder(maxorder);
1201 
1202  data.p = this->MaxOrder();
1203  data.HSize = 2.*this->InnerRadius();
1204  error.Fill(0.);
1205  int npoints = intrule->NPoints(), ip;
1206  for(ip=0;ip<npoints;ip++){
1207  intrule->Point(ip,intpoint,weight);
1208  this->ComputeShape(intpoint, data.x, data.jacobian, data.axes, data.detjac, data.jacinv, data.phi, data.phi, data.dphix);
1209  weight *= fabs(data.detjac);
1210  this->ComputeSolution(intpoint, data.phi, data.dphix, data.axes, data.sol, data.dsol);
1211  material->ContributeErrors(data,weight,error,errorid);
1212  }
1213  intrule->SetOrder(prevorder);
1214 }
1215 
1216 TPZVec<STATE> TPZInterpolationSpace::IntegrateSolution(int variable) const {
1217  TPZMaterial * material = Material();
1218  TPZVec<STATE> result;
1219  if(!material){
1220  PZError << "Error at " << __PRETTY_FUNCTION__ << " : no material for this element\n";
1221  return result;
1222  }
1223  if (!this->Reference()){
1224  PZError << "Error at " << __PRETTY_FUNCTION__ << " : no reference element\n";
1225  return result;
1226  }
1227  const int dim = this->Dimension();
1228  int meshdim = Mesh()->Dimension();
1229  REAL weight;
1230  TPZMaterialData data;
1231  TPZInterpolationSpace *thisnonconst = (TPZInterpolationSpace *) this;
1232 
1233  TPZInterpolationSpace *effective = thisnonconst;
1234  TPZTransform<REAL> tr(dim);
1235  if (dim != Mesh()->Dimension()) {
1236  TPZGeoElSide gelside(thisnonconst->Reference(),this->Reference()->NSides()-1);
1237  TPZGeoElSide neighbour = gelside.Neighbour();
1238  while(neighbour != gelside)
1239  {
1240  if (neighbour.Element()->Reference() && neighbour.Element()->Dimension() == Mesh()->Dimension()) {
1241  break;
1242  }
1243  neighbour = neighbour.Neighbour();
1244  }
1245  if (neighbour == gelside) {
1246  DebugStop();
1247  }
1248  gelside.SideTransform3(neighbour, tr);
1249  TPZTransform<REAL> tr2 = neighbour.Element()->SideToSideTransform(neighbour.Side(), neighbour.Element()->NSides()-1);
1250  tr = tr2.Multiply(tr);
1251  effective = dynamic_cast<TPZInterpolationSpace *> (neighbour.Element()->Reference());
1252  material = effective->Material();
1253  }
1254 
1255  if(!effective) DebugStop();
1256  TPZMaterialData data2d;
1257  thisnonconst->InitMaterialData(data2d);
1258  effective->InitMaterialData(data);
1259  data.fNeedsSol = true;
1260 
1261  TPZManVector<REAL, 3> intpoint(dim,0.);
1262  const int varsize = material->NSolutionVariables(variable);
1263  TPZManVector<STATE,3> value(varsize,0.);
1264 
1265  const TPZIntPoints &intrule = this->GetIntegrationRule();
1266  int npoints = intrule.NPoints(), ip, iv;
1267  TPZManVector<STATE> sol(varsize);
1268  for(ip=0;ip<npoints;ip++) {
1269  intrule.Point(ip,intpoint,weight);
1270 
1271  TPZManVector<REAL,3> intpointTR(meshdim,0.);
1272  tr.Apply(intpoint, intpointTR);
1273  //Tiago: Next call is performed only for computing detcaj. The previous method (Solution) has already computed jacobian.
1274  // It means that the next call would not be necessary if I wrote the whole code here.
1275  thisnonconst->Reference()->Jacobian(intpoint, data2d.jacobian, data2d.axes, data2d.detjac, data2d.jacinv);
1276  effective->Reference()->Jacobian(intpointTR, data.jacobian, data.axes, data.detjac, data.jacinv);
1277  data.intLocPtIndex = ip;
1278 
1279  effective->ComputeRequiredData(data, intpointTR);
1280 
1281  sol.Fill(0.);
1282  material->Solution(data, variable, sol);
1283  weight *= fabs(data2d.detjac);
1284  for(iv = 0; iv < varsize; iv++) {
1285 #ifdef STATE_COMPLEX
1286  value[iv] += sol[iv].real()*weight;
1287 #else
1288 
1289  value[iv] += sol[iv]*weight;
1290 #endif
1291  }//for iv
1292  }//for ip
1293  return value;
1294 }//method
1295 
1296 
1297 void TPZInterpolationSpace::Integrate(int variable, TPZVec<STATE> & value)//AQUIFRAN
1298 {
1299  value = IntegrateSolution(variable);
1300 }
1301 
1302 //void TPZInterpolationSpace::IntegrateSolution(TPZVec<STATE> & value){
1303 // TPZMaterial * material = Material();
1304 // if(!material){
1305 // PZError << "Error at " << __PRETTY_FUNCTION__ << " : no material for this element\n";
1306 // return;
1307 // }
1308 // if (!this->Reference()){
1309 // PZError << "Error at " << __PRETTY_FUNCTION__ << " : no reference element\n";
1310 // return;
1311 // }
1312 // const int dim = this->Dimension();
1313 // REAL weight;
1314 // TPZMaterialData data;
1315 // this->InitMaterialData(data);
1316 //
1317 // TPZManVector<REAL, 3> intpoint(dim,0.);
1318 // const int varsize = material->NStateVariables();
1319 // value.Resize(varsize);
1320 // value.Fill(0.);
1321 //
1322 // const TPZIntPoints &intrule = this->GetIntegrationRule();
1323 // int npoints = intrule.NPoints(), ip, iv;
1324 //
1325 // for(ip=0;ip<npoints;ip++){
1326 // intrule.Point(ip,intpoint,weight);
1327 // data.sol.Fill(0.);
1328 // this->ComputeSolution(intpoint, data.sol, data.dsol, data.axes);
1329 // //Tiago: Next call is performet only for computing detcaj. The previous method (Solution) has already computed jacobian.
1330 // // It means that the next call would not be necessary if I write the whole code here.
1331 // this->Reference()->Jacobian(intpoint, data.jacobian, data.axes, data.detjac, data.jacinv);
1332 // weight *= fabs(data.detjac);
1333 // for(iv = 0; iv < varsize; iv++){
1334 // value[iv] += data.sol[0][iv]*(STATE)weight;
1335 // }//for iv
1336 // }//for ip
1337 //}//method
1338 
1339 void TPZInterpolationSpace::ProjectFlux(TPZElementMatrix &ek, TPZElementMatrix &ef) {
1340 
1341  TPZMaterial * material = Material();
1342  if(!material){
1343  std::stringstream sout;
1344  sout << "Exiting ProjectFlux: no material for this element\n";
1345  Print(sout);
1346  LOGPZ_ERROR(logger,sout.str());
1347  ek.Reset();
1348  ef.Reset();
1349  return;
1350  }
1351 
1352  int num_flux = material->NFluxes();
1353  int dim = Dimension();
1354  int nshape = NShapeF();
1355  int ncon = NConnects();
1356  const TPZIntPoints &intrule = GetIntegrationRule();
1357 
1358  int numeq = nshape;
1359  ek.fMat.Resize(numeq,numeq);
1360  ek.fBlock.SetNBlocks(ncon);
1361  ef.fMat.Resize(numeq,num_flux);
1362  ef.fBlock.SetNBlocks(ncon);
1363 
1364  for(int i=0; i<ncon; ++i){
1365  (ef.fConnect)[i] = ConnectIndex(i);
1366  (ek.fConnect)[i] = ConnectIndex(i);
1367  }
1368 
1369  TPZMaterialData data;
1370  this->InitMaterialData(data);
1371 
1372  //TPZManVector<REAL> flux(num_flux,1);
1373  TPZManVector<STATE> flux(num_flux,1);
1374  TPZManVector<REAL,3> intpoint(dim);
1375  REAL weight = 0.;
1376  for(int int_ind = 0; int_ind < intrule.NPoints(); ++int_ind){
1377 
1378  intrule.Point(int_ind,intpoint,weight);
1379  this->ComputeShape(intpoint, data.x, data.jacobian, data.axes, data.detjac, data.jacinv, data.phi, data.dphi, data.dphix);
1380  weight *= fabs(data.detjac);
1381  this->ComputeSolution(intpoint, data.phi, data.dphix, data.axes, data.sol, data.dsol);
1382 
1383  material->Flux(data.x,data.sol[0],data.dsol[0],data.axes,flux);
1384  for(int in=0; in<nshape; in++){
1385  for(int ifl=0; ifl<num_flux; ifl++){
1386  (ef.fMat)(in,ifl) += flux[ifl]*(STATE)data.phi(in,0)*(STATE)weight;
1387  }//for ifl
1388  for(int jn = 0; jn<nshape; jn++){
1389  (ek.fMat)(in,jn) += data.phi(in,0)*data.phi(jn,0)*weight;
1390  }//for jn
1391  }//for in
1392  }//for int_ind
1393 }//method
1394 
1396 void TPZInterpolationSpace::Write(TPZStream &buf, int withclassid) const
1397 {
1398  TPZCompEl::Write(buf,withclassid);
1399  buf.Write(&fPreferredOrder,1);
1400 }
1401 
1402 
1403 void TPZInterpolationSpace::MinMaxSolutionValues(TPZVec<STATE> &min, TPZVec<STATE> &max){
1404 #ifndef STATE_COMPLEX
1405 
1406  const int dim = Dimension();
1407  TPZManVector<REAL,3> intpoint(dim,0.);
1408 
1409  TPZAutoPointer<TPZIntPoints> intrule = GetIntegrationRule().Clone();
1410  TPZManVector<int,3> prevorder(dim,0);
1411  intrule->GetOrder(prevorder);
1412 
1413  TPZManVector<int,3> maxorder(dim,intrule->GetMaxOrder());
1414  intrule->SetOrder(maxorder);
1415 
1416  TPZSolVec sol;
1417  TPZGradSolVec dsol;
1418  TPZFNMatrix<9> axes(3,3,0.);
1419  REAL weight;
1420 
1421  int intrulepoints = intrule->NPoints();
1422  intrule->Point(0,intpoint,weight);
1423  this->ComputeSolution(intpoint, sol, dsol, axes);
1424  min = sol[0];
1425  max = sol[0];
1426  const int nvars = sol.NElements();
1427  for(int int_ind = 1; int_ind < intrulepoints; int_ind++){
1428  intrule->Point(int_ind,intpoint,weight);
1429  this->ComputeSolution(intpoint, sol, dsol, axes);
1430  for(int iv = 0; iv < nvars; iv++){
1431  if (sol[0][iv] < min[iv]) min[iv] = sol[0][iv];
1432  if (sol[0][iv] > max[iv]) max[iv] = sol[0][iv];
1433  }//iv
1434  }//loop over integratin points
1435 
1436  intrule->SetOrder(prevorder);
1437 #else
1438  DebugStop();
1439 #endif
1440 
1441 }//void
1442 
1443 
1444 void TPZInterpolationSpace::BuildTransferMatrix(TPZInterpolationSpace &coarsel, TPZTransform<> &t, TPZTransfer<STATE> &transfer){
1445  // accumulates the transfer coefficients between the current element and the
1446  // coarse element into the transfer matrix, using the transformation t
1447  TPZGeoEl *ref = Reference();
1448  int locnod = NConnects();
1449  int cornod = coarsel.NConnects();
1450  int locmatsize = NShapeF();
1451  int cormatsize = coarsel.NShapeF();
1452 
1453  // compare interpolation orders
1454  // the minimum interpolation order of this needs to be larger than the maximum interpolation order of coarse
1455 
1456  int mymaxorder = MaxOrder();
1457  int ic;
1458  int coarsemaxorder = this->MaxOrder();
1459  if(coarsemaxorder > mymaxorder) {
1460  std::stringstream sout;
1461  sout << "Exiting BuildTransferMatrix - compute the transfer matrix coarse "
1462  << coarsemaxorder << " me " << mymaxorder << std::endl;
1463  LOGPZ_ERROR(logger,sout.str());
1464  return;
1465  }
1466  TPZStack<int64_t> connectlistcoarse;
1467  TPZStack<int> corblocksize;
1468  TPZStack<int> dependencyordercoarse;
1469  connectlistcoarse.Resize(0);
1470  dependencyordercoarse.Resize(0);
1471  corblocksize.Resize(0);
1472  for(ic=0; ic<coarsel.NConnects(); ic++) connectlistcoarse.Push(coarsel.ConnectIndex(ic));
1473  coarsel.BuildConnectList(connectlistcoarse);
1474  TPZConnect::BuildDependencyOrder(connectlistcoarse,dependencyordercoarse,*coarsel.Mesh());
1475 
1476  // cornod = number of connects associated with the coarse element
1477  cornod = connectlistcoarse.NElements();
1478  int nvar = coarsel.Material()->NStateVariables();
1479 
1480  // number of blocks is cornod
1481  TPZBlock<REAL> corblock(0,cornod);
1482  int in;
1483 
1484  cormatsize = 0;
1485  for(in=0;in<cornod; in++) {
1486  int c = connectlistcoarse[in];
1487  unsigned int blsize = coarsel.Mesh()->ConnectVec()[c].NDof(*(coarsel.Mesh()))/nvar;
1488 #ifdef PZDEBUG
1489  TPZConnect &con = coarsel.Mesh()->ConnectVec()[c];
1490  if(con.NShape() != blsize)
1491  {
1492  DebugStop();
1493  }
1494 #endif
1495  corblock.Set(in,blsize);
1496  corblocksize.Push(blsize);
1497  cormatsize += blsize;
1498  }
1499  corblock.Resequence();
1500 
1501  // REAL loclocmatstore[500] = {0.};
1502  // loclocmat is the inner product of the shape functions of the local element
1503  // loccormat is the inner product of the shape functions with the shape functions
1504  // of the coarse element, both dependent and independent
1505  TPZFNMatrix<500,STATE> loclocmat(locmatsize,locmatsize);
1506  TPZFNMatrix<500,STATE> loccormat(locmatsize,cormatsize);
1507  loclocmat.Zero();
1508  loccormat.Zero();
1509 
1510  TPZAutoPointer<TPZIntPoints> intrule = GetIntegrationRule().Clone();
1511  int dimension = Dimension();
1512 
1513  TPZManVector<int> prevorder(dimension),order(dimension);
1514  intrule->GetOrder(prevorder);
1515 
1516 
1517  // compute the interpolation order of the shapefunctions squared
1518  int dim;
1519  for(dim=0; dim<dimension; dim++) {
1520  order[dim] = mymaxorder*2;
1521  }
1522  intrule->SetOrder(order);
1523 
1524  TPZBlock<REAL> locblock(0,locnod);
1525 
1526  for(in = 0; in < locnod; in++) {
1527  TPZConnect &c = Connect(in);
1528  unsigned int nshape = c.NShape();
1529 #ifdef PZDEBUG
1530  if(NConnectShapeF(in, c.Order()) != nshape)
1531  {
1532  DebugStop();
1533  }
1534 #endif
1535  locblock.Set(in,nshape);
1536  }
1537  locblock.Resequence();
1538 
1539  REAL locphistore[50]={0.},locdphistore[150]={0.};
1540  TPZFMatrix<REAL> locphi(locmatsize,1,locphistore,50);
1541  TPZFMatrix<REAL> locdphi(dimension,locmatsize,locdphistore,150);
1542  locphi.Zero();
1543  locdphi.Zero();
1544  // derivative of the shape function
1545  // in the master domain
1546 
1547  TPZFMatrix<REAL> corphi(cormatsize,1);
1548  TPZFMatrix<REAL> cordphi(dimension,cormatsize);
1549  // derivative of the shape function
1550  // in the master domain
1551 
1552  REAL jacobianstore[9],
1553  axesstore[9];
1554  TPZManVector<REAL> int_point(dimension),
1555  coarse_int_point(dimension);
1556  TPZFMatrix<REAL> jacobian(dimension,dimension,jacobianstore,9),jacinv(dimension,dimension);
1557  TPZFMatrix<REAL> axes(dimension,3,axesstore,9);
1558  TPZManVector<REAL> x(3);
1559  int_point.Fill(0.,0);
1560  REAL jac_det = 1.;
1561  ref->Jacobian( int_point, jacobian , axes, jac_det, jacinv);
1562  REAL multiplier = 1./jac_det;
1563 
1564  int numintpoints = intrule->NPoints();
1565  REAL weight;
1566  int lin,ljn,cjn;
1567 
1568 #ifdef LOG4CXX
1569  if (logger->isDebugEnabled() && coarsel.HasDependency()) {
1570  std::stringstream sout;
1571  coarsel.Print(sout);
1572  int nc = coarsel.NConnects();
1573  for (int ic=0; ic<nc ; ic++) {
1574  coarsel.Connect(ic).Print(*(coarsel.Mesh()),sout);
1575  }
1576  sout << "Connect list coarse " << connectlistcoarse << std::endl;
1577  sout << "dependency order " << dependencyordercoarse << std::endl;
1578  LOGPZ_DEBUG(logger, sout.str())
1579  }
1580 #endif
1581 
1582  for(int int_ind = 0; int_ind < numintpoints; ++int_ind) {
1583 
1584  intrule->Point(int_ind,int_point,weight);
1585  ref->Jacobian( int_point, jacobian , axes, jac_det, jacinv);
1586  ref->X(int_point, x);
1587  Shape(int_point,locphi,locdphi);
1588  weight *= jac_det;
1589  t.Apply(int_point,coarse_int_point);
1590  corphi.Zero();
1591  cordphi.Zero();
1592  coarsel.Shape(coarse_int_point,corphi,cordphi);
1593 
1594 #ifdef LOG4CXX
1595  if (logger->isDebugEnabled() && coarsel.HasDependency()) {
1596  std::stringstream sout;
1597  corphi.Print("Coarse shape functions before expandShapeFunctions",sout);
1598  LOGPZ_DEBUG(logger, sout.str())
1599  }
1600 #endif
1601  coarsel.ExpandShapeFunctions(connectlistcoarse,dependencyordercoarse,corblocksize,corphi,cordphi);
1602 #ifdef LOG4CXX
1603  if (logger->isDebugEnabled() && coarsel.HasDependency()) {
1604  std::stringstream sout;
1605  corphi.Print("Coarse shape functions after expandShapeFunctions",sout);
1606  LOGPZ_DEBUG(logger, sout.str())
1607  }
1608 #endif
1609 
1610  for(lin=0; lin<locmatsize; lin++) {
1611  for(ljn=0; ljn<locmatsize; ljn++) {
1612  loclocmat(lin,ljn) += weight*locphi(lin,0)*locphi(ljn,0)*multiplier;
1613  }
1614  for(cjn=0; cjn<cormatsize; cjn++) {
1615  loccormat(lin,cjn) += weight*locphi(lin,0)*corphi(cjn,0)*multiplier;
1616  }
1617  }
1618  jacobian.Zero();
1619  }
1620  loclocmat.SolveDirect(loccormat,ELDLt);
1621 
1622 #ifdef LOG4CXX
1623  {
1624  std::stringstream sout;
1625  loccormat.Print("Element transfer matrix",sout);
1626  LOGPZ_DEBUG(logger, sout.str())
1627  }
1628 #endif
1629 
1630  for(in=0; in<locnod; in++) {
1631  // int cind = connectlistcoarse[in];
1632  if(Connect(in).HasDependency()) continue;
1633  int64_t locblocknumber = Connect(in).SequenceNumber();
1634  int locblocksize = locblock.Size(in);
1635  int64_t locblockpos = locblock.Position(in);
1636  TPZStack<int> locblockvec;
1637  TPZStack<int> globblockvec;
1638  int64_t numnonzero = 0, jn;
1639  // if(transfer.HasRowDefinition(locblocknumber)) continue;
1640 
1641  for(jn = 0; jn<cornod; jn++) {
1642  int corblocksize = corblock.Size(jn);
1643  int64_t corblockpos = corblock.Position(jn);
1644  int64_t cind = connectlistcoarse[jn];
1645  TPZConnect &con = coarsel.Mesh()->ConnectVec()[cind];
1646  if(con.HasDependency()) continue;
1647  int64_t corblocknumber = con.SequenceNumber();
1648  if(locblocksize == 0 || corblocksize == 0) continue;
1649  TPZFMatrix<STATE> smalll(locblocksize,corblocksize,0.);
1650  loccormat.GetSub(locblockpos,corblockpos, locblocksize,corblocksize,smalll);
1651  REAL tol = Norm(smalll);
1652  if(tol >= 1.e-10) {
1653  locblockvec.Push(jn);
1654  globblockvec.Push(corblocknumber);
1655  numnonzero++;
1656  }
1657  }
1658  if(transfer.HasRowDefinition(locblocknumber)) continue;
1659  transfer.AddBlockNumbers(locblocknumber,globblockvec);
1660  int64_t jnn;
1661  for(jnn = 0; jnn<numnonzero; jnn++) {
1662  jn = locblockvec[jnn];
1663  int corblocksize = corblock.Size(jn);
1664  int64_t corblockpos = corblock.Position(jn);
1665  if(corblocksize == 0 || locblocksize == 0) continue;
1666  TPZFMatrix<STATE> smalll(locblocksize,corblocksize,0.);
1667  loccormat.GetSub(locblockpos,corblockpos,locblocksize,corblocksize,smalll);
1668  transfer.SetBlockMatrix(locblocknumber,globblockvec[jnn],smalll);
1669  }
1670  }
1671  intrule->SetOrder(prevorder);
1672 }
1673 
1674 
1675 void TPZInterpolationSpace::ExpandShapeFunctions(TPZVec<int64_t> &connectlist, TPZVec<int> &dependencyorder, TPZVec<int> &blocksizes, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphix) {
1676  int64_t numblocks = connectlist.NElements();
1677  TPZCompMesh &mesh = *Mesh();
1678  int nhandled=0;
1679  int current_order = 0;
1680  int current_block =0;
1681  while(nhandled < numblocks) {
1682  if(dependencyorder[current_block] == current_order) {
1683  nhandled++;
1684  int64_t cind = connectlist[current_block];
1685  TPZConnect &con = mesh.ConnectVec()[cind];
1686  con.ExpandShape(cind,connectlist,blocksizes,phi,dphix);
1687  }
1688  current_block++;
1689  if(current_block == numblocks) {
1690  current_block = 0;
1691  current_order++;
1692  }
1693  }
1694 }
1695 
1696 int TPZInterpolationSpace::GetSideOrient(int side){
1697 
1698  TPZGeoEl *gel = this->Reference();
1699  int nsides = gel->NSides();
1700  int nnos = gel->NCornerNodes();
1701  // para elemento quadratico, ie, quadrilatero de 8 nos
1702  if (gel->Type()==EQuadrilateral && nnos>side) {
1703  int sideorient = gel->NormalOrientation(side);
1704  return sideorient;
1705  }
1706 
1707  if(side < nnos || side ==nsides-1) DebugStop();
1708 
1709  int sideorient = gel->NormalOrientation(side);
1710  return sideorient;
1711 }
1712 
1718 void TPZInterpolationSpace::SetSideOrient(int side, int sideorient){
1719  DebugStop();
1720 }
1721 
1722 
1724 void TPZInterpolationSpace::Read(TPZStream &buf, void *context)
1725 {
1726  TPZCompEl::Read(buf,context);
1727  buf.Read(&fPreferredOrder,1);
1728 }
1729 
1731 void TPZInterpolationSpace::Convert2Axes(const TPZFMatrix<REAL> &dphi, const TPZFMatrix<REAL> &jacinv, TPZFMatrix<REAL> &dphidx)
1732 {
1733  int nshape = dphi.Cols();
1734  int dim = dphi.Rows();
1735  dphidx.Resize(dim,nshape);
1736  int ieq;
1737  switch(dim){
1738  case 0:
1739  {
1740 
1741  }
1742  break;
1743  case 1:
1744  {
1745  dphidx = dphi;
1746  dphidx *= jacinv.GetVal(0,0);
1747  }
1748  break;
1749  case 2:
1750  {
1751  for(ieq = 0; ieq < nshape; ieq++) {
1752  dphidx(0,ieq) = jacinv.GetVal(0,0)*dphi.GetVal(0,ieq) + jacinv.GetVal(1,0)*dphi.GetVal(1,ieq);
1753  dphidx(1,ieq) = jacinv.GetVal(0,1)*dphi.GetVal(0,ieq) + jacinv.GetVal(1,1)*dphi.GetVal(1,ieq);
1754  }
1755  }
1756  break;
1757  case 3:
1758  {
1759  for(ieq = 0; ieq < nshape; ieq++) {
1760  dphidx(0,ieq) = jacinv.GetVal(0,0)*dphi.GetVal(0,ieq) + jacinv.GetVal(1,0)*dphi.GetVal(1,ieq) + jacinv.GetVal(2,0)*dphi.GetVal(2,ieq);
1761  dphidx(1,ieq) = jacinv.GetVal(0,1)*dphi.GetVal(0,ieq) + jacinv.GetVal(1,1)*dphi.GetVal(1,ieq) + jacinv.GetVal(2,1)*dphi.GetVal(2,ieq);
1762  dphidx(2,ieq) = jacinv.GetVal(0,2)*dphi.GetVal(0,ieq) + jacinv.GetVal(1,2)*dphi.GetVal(1,ieq) + jacinv.GetVal(2,2)*dphi.GetVal(2,ieq);
1763  }
1764  }
1765  break;
1766  default:
1767  {
1768  PZError << "Error at " << __PRETTY_FUNCTION__ << " please implement the " << dim << "d Jacobian and inverse\n";
1769  }
1770  } //switch
1771 
1772 }
1773 
1774 int TPZInterpolationSpace::ClassId() const{
1775  return Hash("TPZInterpolationSpace") ^ TPZCompEl::ClassId() << 1;
1776 }
TPZInterpolationSpace::~TPZInterpolationSpace
virtual ~TPZInterpolationSpace()
Default destructor.
Definition: pzinterpolationspace.cpp:52
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZCompMesh::ElementSolution
TPZFMatrix< STATE > & ElementSolution()
Access method for the element solution vectors.
Definition: pzcmesh.h:225
TPZIntPoints::GetOrder
virtual void GetOrder(TPZVec< int > &ord) const =0
Gets the order of the integration rule for each dimension of the master element.
TPZCompMesh::Reference
TPZGeoMesh * Reference() const
Returns a pointer to the geometrical mesh associated.
Definition: pzcmesh.h:209
EQuadrilateral
Definition: pzeltype.h:57
TPZInterpolationSpace::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: pzinterpolationspace.cpp:1774
fabs
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
TPZInterpolationSpace::CreateInterfaces
void CreateInterfaces(bool BetweenContinuous=false)
Create interfaces between this and its neighbours.
Definition: pzinterpolationspace.cpp:658
pzquad.h
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
TPZInterpolationSpace::GetSideOrient
virtual int GetSideOrient(int side)
It returns the normal orientation of the reference element by the side. Only side that has dimension ...
Definition: pzinterpolationspace.cpp:1696
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZBlock::Position
int Position(const int block_diagonal) const
Returns the position of first element block dependent on matrix diagonal.
Definition: pzblock.h:177
TPZIntPoints::NPoints
virtual int NPoints() const =0
Returns number of points for the cubature rule related.
TPZCompElSide::LowerLevelElementList
TPZCompElSide LowerLevelElementList(int onlyinterpolated)
Returns all connected elements which have level lower to the current element.
Definition: pzcompel.cpp:824
TPZCompElSide
Implements computational element and a side. Computational Element.
Definition: pzcompel.h:632
TPZBlock::Set
int Set(const int index, const int dim, const int pos=-1)
Modifies existing block dimensions or creates a new block with given index.
Definition: pzblock.cpp:104
TPZConnect
Represents a set of shape functions associated with a computational element/side. Computational Eleme...
Definition: pzconnect.h:30
tpzchangeel.h
Contains the TPZChangeEl class. It is a special map.
ELDLt
Definition: pzmatrix.h:52
TPZMaterial::Errors
virtual void Errors(TPZMaterialData &data, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &errors)
Definition: TPZMaterial.h:496
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZCompEl::Type
virtual MElementType Type()
Return the type of the element.
Definition: pzcompel.cpp:195
TPZInterpolationSpace::InitMaterialData
virtual void InitMaterialData(TPZMaterialData &data)
Initialize a material data and its attributes based on element dimension, number of state variables a...
Definition: pzinterpolationspace.cpp:127
TPZGeoEl::MaterialId
int MaterialId() const
Returns the material index of the element.
Definition: pzgeoel.h:250
TPZInterpolationSpace::ExistsInterface
int ExistsInterface(TPZGeoElSide geosd)
Verify existence of interface.
Definition: pzinterpolationspace.cpp:964
TPZTransfer::SetBlockMatrix
void SetBlockMatrix(int row, int col, TPZFMatrix< TVar > &mat)
Sets the row,col block equal to matrix mat if row col was not specified by AddBlockNumbers, an error will be issued and exit.
Definition: pztransfer.cpp:153
TPZCompEl::BuildConnectList
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
ETriangle
Definition: pzeltype.h:56
TPZVec::resize
virtual void resize(const int64_t newsize)
Definition: pzvec.h:213
TPZInterpolationSpace::NConnectShapeF
virtual int NConnectShapeF(int icon, int order) const =0
Returns the number of shapefunctions associated with a connect.
TPZMaterial::ErrorsHdiv
virtual void ErrorsHdiv(TPZMaterialData &data, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values)
Definition: TPZMaterial.h:518
dimension
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
TPZElementMatrix::fConnect
TPZStack< int64_t > fConnect
Vector of pointers to TPZConnect objects.
Definition: pzelmat.h:39
TPZInterpolationSpace::Integrate
virtual void Integrate(int variable, TPZVec< STATE > &value) override
Integrates a variable over the element.
Definition: pzinterpolationspace.cpp:1297
TPZTransfer::AddBlockNumbers
void AddBlockNumbers(int row, TPZVec< int > &colnumbers)
Will specify the sparsity pattern of row.
Definition: pztransfer.cpp:121
TPZMaterialData::dphi
TPZFNMatrix< 660, REAL > dphi
values of the derivative of the shape functions over the master element
Definition: pzmaterialdata.h:57
TPZCompElSide::Reference
TPZGeoElSide Reference() const
Reference to the geometric element.
Definition: pzcompel.cpp:748
TPZElementMatrix::fMat
TPZFNMatrix< 1000, STATE > fMat
Pointer to a blocked matrix object.
Definition: pzelmat.h:41
TPZInterpolationSpace::ComputeIntegrationOrder
virtual int ComputeIntegrationOrder() const override
Compute integration order according to ... .
Definition: pzinterpolationspace.cpp:85
TPZElementMatrix::fOneRestraints
std::list< TPZOneShapeRestraint > fOneRestraints
list of one degree of freedom restraints
Definition: pzelmat.h:56
TPZGeoEl::BuildTransform2
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
TPZCompEl::Dimension
virtual int Dimension() const =0
Dimension of the element.
TPZGeoEl::NCornerNodes
virtual int NCornerNodes() const =0
Returns the number of corner nodes of the element.
pzinterpolationspace.h
Contains declaration of TPZInterpolationSpace class which implements the interface for interpolated c...
TPZMaterial::Dimension
virtual int Dimension() const =0
Returns the integrable dimension of the material.
TPZMaterial::NStateVariables
virtual int NStateVariables() const =0
Returns the number of state variables associated with the material.
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZMaterialData::fShapeType
MShapeFunctionType fShapeType
Definition: pzmaterialdata.h:43
TPZMaterialData::jacinv
TPZFNMatrix< 9, REAL > jacinv
value of the inverse of the jacobian at the integration point
Definition: pzmaterialdata.h:67
TPZGeoElSide
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
TPZInterpolationSpace::CalcResidual
virtual void CalcResidual(TPZElementMatrix &ef) override
Only computes the element residual.
Definition: pzinterpolationspace.cpp:377
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
LOGPZ_WARN
#define LOGPZ_WARN(A, B)
Define log for warnings.
Definition: pzlog.h:91
TPZGeoElSide::Neighbour
TPZGeoElSide Neighbour() const
Definition: pzgeoel.h:754
TPZGeoEl::NSides
virtual int NSides() const =0
Returns the number of connectivities of the element.
TPZInterpolationSpace::GetIntegrationRule
virtual const TPZIntPoints & GetIntegrationRule() const override=0
Returns a reference to an integration rule suitable for integrating the interior of the element...
TPZMaterial::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)=0
It computes a contribution to the stiffness matrix and load vector at one integration point...
TPZElementMatrix::fBlock
TPZBlock< STATE > fBlock
Block structure associated with fMat.
Definition: pzelmat.h:43
TPZMaterialData::xParametric
TPZManVector< REAL, 3 > xParametric
value of the coordinate at the integration point corresponding to the x-parametric coordinate (master...
Definition: pzmaterialdata.h:73
TPZGeoEl::SideDimension
virtual int SideDimension(int side) const =0
Return the dimension of side.
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZMaterialData::jacobian
TPZFNMatrix< 9, REAL > jacobian
value of the jacobian at the integration point
Definition: pzmaterialdata.h:65
pztransfer.h
Contains the TPZTransfer class which implements a rectangular sparse block matrix.
TPZInterpolationSpace::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzinterpolationspace.cpp:1724
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZMaterialData::p
int p
maximum polinomial order of the shape functions
Definition: pzmaterialdata.h:75
TPZMatrix::GetSub
virtual int GetSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize, const int64_t colSize, TPZFMatrix< TVar > &Target) const
Gets submatrix storing it on Target.
Definition: pzmatrix.cpp:602
TPZInterpolationSpace::TPZInterpolationSpace
TPZInterpolationSpace()
Default constructor.
Definition: pzinterpolationspace.cpp:22
LOGPZ_FATAL
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Definition: pzlog.h:95
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZCompEl::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: pzcompel.cpp:1129
TPZIntPoints
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
TPZInterpolationSpace::SetIntegrationRule
virtual void SetIntegrationRule(int order) override
Definition: pzinterpolationspace.h:106
error
void error(char *string)
Definition: testShape.cc:7
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
TPZGeoEl::CenterPoint
virtual void CenterPoint(int side, TPZVec< REAL > &masscent) const =0
It returns the coordinates from the center of the side of the element in the element coordinate space...
TPZInterpolationSpace::BuildTransferMatrix
void BuildTransferMatrix(TPZInterpolationSpace &coarsel, TPZTransform<> &t, TPZTransfer< STATE > &transfer)
Accumulates the transfer coefficients between the current element and the coarse element into the t...
Definition: pzinterpolationspace.cpp:1444
TPZInterpolationSpace::Shape
virtual void Shape(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphidxi)=0
Computes the shape function set at the point x.
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZCompEl::Material
virtual TPZMaterial * Material() const
Identify the material object associated with the element.
Definition: pzcompel.cpp:959
TPZCompElSide::HigherLevelElementList
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
TPZInterpolationSpace::ProjectFlux
void ProjectFlux(TPZElementMatrix &ek, TPZElementMatrix &ef) override
Integrate the solution over the element.
Definition: pzinterpolationspace.cpp:1339
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZMaterial::Flux
virtual void Flux(TPZVec< REAL > &x, TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux)
Computes the value of the flux function to be used by ZZ error estimator.
Definition: TPZMaterial.h:255
TPZMaterial::IntegrationRuleOrder
virtual int IntegrationRuleOrder(int elPMaxOrder) const
Gets the order of the integration rule necessary to integrate an element with polinomial order p...
Definition: TPZMaterial.cpp:349
TPZConnect::Order
unsigned char Order() const
Access function to return the order associated with the connect.
Definition: pzconnect.h:217
TPZConnect::SequenceNumber
int64_t SequenceNumber() const
Returns the Sequence number of the connect object.
Definition: pzconnect.h:158
TPZMaterial::NEvalErrors
virtual int NEvalErrors()
Returns the number of norm errors. Default is 3: energy, L2 and H1.
Definition: TPZMaterial.h:524
TPZMaterialData::XCenter
TPZManVector< REAL, 3 > XCenter
value of the coordinate at the center of the element
Definition: pzmaterialdata.h:87
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
TPZCompEl::Read
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzcompel.cpp:736
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
TPZInterpolationSpace::RemoveInterfaces
void RemoveInterfaces()
Remove interfaces connected to this element.
Definition: pzinterpolationspace.cpp:977
TPZCompElSide::EqualLevelElementList
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
LOGPZ_INFO
#define LOGPZ_INFO(A, B)
Define log for informations.
Definition: pzlog.h:89
TPZInterpolationSpace::ComputeShape
virtual void ComputeShape(TPZVec< REAL > &intpoint, TPZVec< REAL > &X, TPZFMatrix< REAL > &jacobian, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< REAL > &dphidx)
Compute shape functions based on master element in the classical FEM manner.
Definition: pzinterpolationspace.cpp:96
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZTransfer::HasRowDefinition
int HasRowDefinition(int row)
Returns 1 if the row is defined (i.e. has column entries)
Definition: pztransfer.cpp:114
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
TPZMaterial::NFluxes
virtual int NFluxes()
Returns the number of components which form the flux function.
Definition: TPZMaterial.h:183
TPZMaterial::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data)
Fill material data parameter with necessary requirements for the.
Definition: TPZMaterial.cpp:81
TPZInterpolationSpace::CreateInterface
TPZInterfaceElement * CreateInterface(int side, bool BetweenContinuous=false)
Create an interface between this and the neighbour by side side.
Definition: pzinterpolationspace.cpp:710
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZInterfaceEl.h
Contains declaration of TPZInterfaceElement class which computes the contribution over an interface b...
TPZIntPoints::SetOrder
virtual void SetOrder(TPZVec< int > &ord, int type=0)=0
Sets the order of the cubature rule.
TPZInterpolationSpace::Convert2Axes
static void Convert2Axes(const TPZFMatrix< REAL > &dphi, const TPZFMatrix< REAL > &jacinv, TPZFMatrix< REAL > &dphidx)
convert a shapefunction derivative in xi-eta to a function derivative in axes
Definition: pzinterpolationspace.cpp:1731
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZCompEl::Mesh
TPZCompMesh * Mesh() const
Return a pointer to the grid of the element.
Definition: pzcompel.cpp:288
ELU
Definition: pzmatrix.h:52
TPZCompMesh::Dimension
int Dimension() const
Returns the dimension of the simulation.
Definition: pzcmesh.h:148
TPZGeoEl::NumInterfaces
int NumInterfaces()
Returns number of TPZInterfaceElement pointing to this.
Definition: pzgeoel.h:94
TPZGeoEl::Type
virtual MElementType Type() const =0
Returns the element type acording to pzeltype.h.
TPZMaterial::ContributeErrors
virtual void ContributeErrors(TPZMaterialData &data, REAL weight, TPZVec< REAL > &nk, int &errorid)
Definition: TPZMaterial.h:538
TPZGeoEl::Jacobian
void Jacobian(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &jac, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Compute a decomposition of the gradient of the mapping function, as a rotation matrix (Jacobian) and ...
Definition: pzgeoel.cpp:1144
TPZCompMesh::FindMaterial
TPZMaterial * FindMaterial(int id)
Find the material with identity id.
Definition: pzcmesh.cpp:297
TPZGeoEl::Reference
TPZCompEl * Reference() const
Return a pointer to the element referenced by the geometric element.
Definition: pzgeoel.h:750
TPZConnect::NShape
unsigned int NShape() const
Definition: pzconnect.h:151
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZConnect::HasDependency
int HasDependency() const
Returns whether exist dependecy information.
Definition: pzconnect.h:292
TPZInterpolationSpace::ComputeError
virtual void ComputeError(int errorid, TPZVec< REAL > &error) override
Computes the element error estimator.
Definition: pzinterpolationspace.cpp:1180
TPZGeoEl::NormalOrientation
int NormalOrientation(int side)
Determine the orientation of the normal vector comparing the ids of the neighbouring elements...
Definition: pzgeoel.cpp:2370
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZGeoEl::CreateBCGeoEl
virtual TPZGeoEl * CreateBCGeoEl(int side, int bc)=0
Method which creates a geometric element on the side of an existing element.
TPZMaterialData::axes
TPZFNMatrix< 9, REAL > axes
axes indicating the directions of the derivatives of the shapefunctions
Definition: pzmaterialdata.h:63
TPZCompMesh::Block
const TPZBlock< STATE > & Block() const
Access the block structure of the solution vector.
Definition: pzcmesh.h:213
TPZElementMatrix::fMesh
TPZCompMesh * fMesh
Definition: pzelmat.h:36
TPZInterpolationSpace::Solution
virtual void Solution(TPZVec< REAL > &qsi, int var, TPZVec< STATE > &sol) override
Post processing method which computes the solution for the var post processed variable.
Definition: pzinterpolationspace.cpp:508
TPZCompElSide::Exists
int Exists() const
Verifies if the object is non null (initialized)
Definition: pzcompel.h:697
TPZCompElDisc.h
Contains declaration of TPZCompelDisc class which implements a computational element for discontinuou...
TPZCompMesh::ConnectVec
TPZAdmChunkVector< TPZConnect > & ConnectVec()
Return a reference to the connect pointers vector.
Definition: pzcmesh.h:198
TPZCompEl::Index
int64_t Index() const
Returns element index of the mesh fELementVec list.
Definition: pzcompel.h:821
pzmaterialdata.h
Contains the TPZMaterialData class which implements an interface between TPZCompEl::CalcStiff and TPZ...
TPZMaterialData::HSize
REAL HSize
measure of the size of the element
Definition: pzmaterialdata.h:83
TPZInterpolationSpace::InitializeElementMatrix
virtual void InitializeElementMatrix(TPZElementMatrix &ek, TPZElementMatrix &ef)
Initialize element matrix in which is computed CalcStiff.
Definition: pzinterpolationspace.cpp:429
TPZMaterialData::intGlobPtIndex
int intGlobPtIndex
global point index
Definition: pzmaterialdata.h:115
TPZGeoEl::Neighbour
virtual TPZGeoElSide Neighbour(int side)=0
Returns a pointer to the neighbour and the neighbourside along side of the current element...
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZBlock::SetNBlocks
int SetNBlocks(const int num_of_blocks)
Sets number of blocks on diagonal matrix.
Definition: pzblock.cpp:91
TPZCompEl::ComputeSolution
virtual void ComputeSolution(TPZVec< REAL > &qsi, TPZMaterialData &data)
Definition: pzcompel.h:462
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZCompEl::Connect
virtual TPZConnect & Connect(int i) const
Returns a pointer to the ith node.
Definition: pzcompel.cpp:298
TPZInterpolationSpace::VectorialProd
void VectorialProd(TPZVec< REAL > &ivec, TPZVec< REAL > &jvec, TPZVec< REAL > &kvec, bool unitary=false)
Computes the vectorial product of two vectors and normalize the result if unitary is set to true...
Definition: pzinterpolationspace.cpp:288
TPZInterpolationSpace::InterpolateSolution
void InterpolateSolution(TPZInterpolationSpace &coarsel)
Interpolates the solution into the degrees of freedom nodes from the degrees of freedom nodes from th...
Definition: pzinterpolationspace.cpp:542
TPZMaterialData::intLocPtIndex
int intLocPtIndex
Index of the current integration point being evaluated.
Definition: pzmaterialdata.h:112
TPZInterpolationSpace::GetPreferredOrder
virtual int GetPreferredOrder()
Returns the prefered order for the element.
Definition: pzinterpolationspace.h:304
TPZConnect::NState
unsigned char NState() const
Number of state variables associated with the connect.
Definition: pzconnect.h:146
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZCompEl::ConnectIndex
virtual int64_t ConnectIndex(int i) const =0
Returns the index of the ith connectivity of the element.
TPZInterfaceElement
Computes the contribution over an interface between two discontinuous elements. Computational Element...
Definition: TPZInterfaceEl.h:30
TPZTransform::Multiply
TPZTransform< T > Multiply(TPZTransform< T > &right)
Multiply the transformation object (to the right) with right (Multiplying matrices) ...
Definition: pztrnsform.cpp:112
TPZElementMatrix
This class associates an element matrix with the coeficients of its contribution in the global stiffn...
Definition: pzelmat.h:30
TPZInterpolationSpace
Implements the interfaces for TPZCompElDisc, TPZInterfaceElement and TPZInterpolatedElement. Computational element.
Definition: pzinterpolationspace.h:17
TPZCompEl::NConnects
virtual int NConnects() const =0
Returns the number of nodes of the element.
TPZGeoEl::IsLinearMapping
virtual bool IsLinearMapping() const
Definition: pzgeoel.h:268
TPZGeoElSide::Element
TPZGeoEl * Element() const
Definition: pzgeoelside.h:162
TPZInterpolationSpace::SetSideOrient
virtual void SetSideOrient(int side, int sideorient)
It set the normal orientation of the element by the side. Only side that has dimension equal to my di...
Definition: pzinterpolationspace.cpp:1718
TPZCompEl::GetShapeRestraints
virtual std::list< TPZOneShapeRestraint > GetShapeRestraints() const
Return a list with the shape restraints.
Definition: pzcompel.h:432
TPZInterpolationSpace::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzinterpolationspace.cpp:1396
TPZInterpolationSpace::IntegrateSolution
virtual TPZVec< STATE > IntegrateSolution(int variable) const override
Integrate a variable over the element.
Definition: pzinterpolationspace.cpp:1216
TPZGeoEl::Dimension
virtual int Dimension() const =0
Returns the dimension of the element.
TPZCompElSide::Element
TPZCompEl * Element() const
Gives a pointer to the reference computational element.
Definition: pzcompel.h:675
TPZGeoEl::Id
int64_t Id() const
Returns the Id of the element.
Definition: pzgeoel.h:223
TPZInterpolationSpace::CalcStiff
virtual void CalcStiff(TPZElementMatrix &ek, TPZElementMatrix &ef) override
Computes the element stiffness matrix and right hand side.
Definition: pzinterpolationspace.cpp:306
TPZCompEl::Reference
TPZGeoEl * Reference() const
Return a pointer to the corresponding geometric element if such exists, return 0 otherwise.
Definition: pzcompel.cpp:1137
TPZInterpolationSpace::NShapeF
virtual int NShapeF() const =0
It returns the shapes number of the element.
TPZInterfaceElement::LeftElement
TPZCompEl * LeftElement() const
Returns the left element from the element interface.
Definition: TPZInterfaceEl.h:188
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZGeoElSide::Reference
TPZCompElSide Reference() const
Returns a pointer to the elementside referenced by the geometric elementside.
Definition: pzgeoelside.cpp:587
TPZGeoEl::X
virtual void X(TPZVec< REAL > &qsi, TPZVec< REAL > &result) const =0
Return the coordinate in real space of the point coordinate in the master element space...
TPZInterpolationSpace::ExpandShapeFunctions
void ExpandShapeFunctions(TPZVec< int64_t > &connectlist, TPZVec< int > &dependencyorder, TPZVec< int > &blocksizes, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
Auxiliary method to expand a vector of shapefunctions and their derivatives to acount for constraints...
Definition: pzinterpolationspace.cpp:1675
TPZIntPoints::GetMaxOrder
virtual int GetMaxOrder() const
Returns the minimum order to integrate polinomials exactly for all implemented cubature rules...
Definition: pzquad.cpp:30
TPZGeoEl::IsGeoBlendEl
virtual bool IsGeoBlendEl() const
Definition: pzgeoel.h:275
TPZGeoElSide::Dimension
int Dimension() const
the dimension associated with the element/side
Definition: pzgeoelside.cpp:593
TPZElementMatrix::Reset
void Reset(TPZCompMesh *mesh=NULL, MType type=Unknown)
Reset the data structure.
Definition: pzelmat.h:59
TPZCompMesh
Implements computational mesh. Computational Mesh.
Definition: pzcmesh.h:47
TPZConnect::BuildDependencyOrder
static void BuildDependencyOrder(TPZVec< int64_t > &connectlist, TPZVec< int > &DependenceOrder, TPZCompMesh &mesh)
This method builds the vector DependenceOrder which indicates in which order constrained nodes need t...
Definition: pzconnect.cpp:556
TPZBlock::Size
int Size(const int block_diagonal) const
Returns block dimension.
Definition: pzblock.h:171
TPZCompEl::Print
virtual void Print(std::ostream &out=std::cout) const
Prints element data.
Definition: pzcompel.cpp:322
TPZCompEl::Solution
virtual void Solution(TPZVec< REAL > &qsi, int var, TPZVec< STATE > &sol)
Calculates the solution - sol - for the variable var at point qsi, where qsi is expressed in terms of...
Definition: pzcompel.cpp:421
TPZMaterial::NumLoadCases
int NumLoadCases()
returns the number of load cases for this material object
Definition: TPZMaterial.h:186
TPZMaterial::Id
int Id() const
Definition: TPZMaterial.h:170
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
TPZCompElSide::Side
int Side() const
Returns the side index.
Definition: pzcompel.h:688
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZCompEl::Write
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzcompel.cpp:726
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZInterpolationSpace::MaxOrder
virtual int MaxOrder()
Returns the max order of interpolation.
Definition: pzinterpolationspace.cpp:55
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZInterpolationSpace::ComputeNormal
virtual void ComputeNormal(TPZMaterialData &data)
Computes the proper normal vector towards the neighbour element.
Definition: pzinterpolationspace.cpp:203
TPZGeoMesh::InterfaceMaterial
int InterfaceMaterial(int leftmaterial, int rightmaterial)
Returns the interface material associated to left and right element materials. If no interface materi...
Definition: pzgmesh.cpp:1548
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZInterfaceElement::RightElement
TPZCompEl * RightElement() const
Returns the right element from the element interface.
Definition: TPZInterfaceEl.h:183
TPZInterpolationSpace::MinMaxSolutionValues
void MinMaxSolutionValues(TPZVec< STATE > &min, TPZVec< STATE > &max)
Returns minimum and maximum values for each state variable.
Definition: pzinterpolationspace.cpp:1403
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
rdt.values
def values
Definition: rdt.py:119
TPZIntPoints::Clone
virtual TPZIntPoints * Clone() const =0
Make a clone of the related cubature rule.
TPZCompElDisc
This class implements a discontinuous element (for use with discontinuous Galerkin). Computational Element.
Definition: TPZCompElDisc.h:35
TPZInterpolationSpace::ComputeRequiredData
virtual void ComputeRequiredData(TPZMaterialData &data, TPZVec< REAL > &qsi)
Compute and fill data with requested attributes.
Definition: pzinterpolationspace.cpp:158
TPZGeoEl::ElementRadius
virtual REAL ElementRadius()
Definition: pzgeoel.cpp:949
TPZInterpolationSpace::InnerRadius
virtual REAL InnerRadius()
Returns the inner radius value.
Definition: pzinterpolationspace.cpp:119
TPZTransfer
Implements rectangular matrix which extends a solution vector of the coarse mesh to a solution vector...
Definition: pzcmesh.h:33
TPZGeoEl::Print
virtual void Print(std::ostream &out=std::cout)
Print all relevant data of the element to cout.
Definition: pzgeoel.cpp:238
TPZInterpolationSpace::fPreferredOrder
int fPreferredOrder
Preferred polynomial order.
Definition: pzinterpolationspace.h:296
TPZMaterialData::fVecShapeIndex
TPZManVector< std::pair< int, int64_t > > fVecShapeIndex
correspondence between normal vector index and index of the shape functions
Definition: pzmaterialdata.h:99
TPZConnect::Print
void Print(const TPZCompMesh &mesh, std::ostream &out=std::cout)
Print the information for the connect element.
Definition: pzconnect.cpp:67
TPZInterpolationSpace::AdjustIntegrationRule
virtual void AdjustIntegrationRule()
Adjust the integration rule according to the polynomial order of shape functions. ...
Definition: pzinterpolationspace.cpp:71
TPZCompEl
Defines the interface of a computational element. Computational Element.
Definition: pzcompel.h:59
EInterface
Definition: pzeltype.h:63
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZCompMesh::GetDefaultOrder
int GetDefaultOrder()
Definition: pzcmesh.h:398
TPZCompEl::HasDependency
virtual int HasDependency()
Returns 1 if the element has at least one dependent node. Returns 0 otherwise.
Definition: pzcompel.cpp:565
TPZFMatrix::GetVal
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
TPZConnect::ExpandShape
void ExpandShape(int64_t cind, TPZVec< int64_t > &connectlist, TPZVec< int > &blocksize, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
Definition: pzconnect.cpp:353
TPZMaterialData::detjac
REAL detjac
determinant of the jacobian
Definition: pzmaterialdata.h:85
TPZTransform::Apply
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZIntPoints::Point
virtual void Point(int i, TPZVec< REAL > &pos, REAL &w) const =0
Returns i-th point at master element and related weight.
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZInterpolationSpace::RemoveInterface
void RemoveInterface(int side)
Remove interface which is neighbour from side side.
Definition: pzinterpolationspace.cpp:1032
TPZBlock::Resequence
int Resequence(const int start=0)
Resequences blocks positioning.
Definition: pzblock.cpp:150
TPZCompEl::fMesh
TPZCompMesh * fMesh
Computational mesh to which the element belongs.
Definition: pzcompel.h:64
TPZInterpolationSpace::Print
virtual void Print(std::ostream &out=std::cout) const override
Prints the relevant data of the element to the output stream.
Definition: pzinterpolationspace.cpp:90
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZInterpolationSpace::EvaluateError
virtual void EvaluateError(std::function< void(const TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv)> func, TPZVec< REAL > &errors, bool store_error) override
Performs an error estimate on the elemen.
Definition: pzinterpolationspace.cpp:1067
Generated on Sun Aug 16 2020 11:17:07 for NeoPZ by   doxygen 1.8.13