diff --git a/ErrorEstimation/CMakeLists.txt b/ErrorEstimation/CMakeLists.txt index 85e6f2e8..db6d1380 100644 --- a/ErrorEstimation/CMakeLists.txt +++ b/ErrorEstimation/CMakeLists.txt @@ -4,8 +4,8 @@ add_library(ErrorEstimationLib ProblemConfig.h TPZHDivErrorEstimator.cpp TPZHDivErrorEstimator.h - TPZHybridH1ErrorEstimator.cpp - TPZHybridH1ErrorEstimator.h + #TPZHybridH1ErrorEstimator.cpp + #TPZHybridH1ErrorEstimator.h TPZMHMHDivErrorEstimator.cpp TPZMHMHDivErrorEstimator.h # Related tools diff --git a/ErrorEstimation/Material/CMakeLists.txt b/ErrorEstimation/Material/CMakeLists.txt index 4a146fa6..0613cc68 100644 --- a/ErrorEstimation/Material/CMakeLists.txt +++ b/ErrorEstimation/Material/CMakeLists.txt @@ -4,14 +4,14 @@ target_sources(ErrorEstimationLib PUBLIC TPZHDivErrorEstimateMaterial.h TPZMixedHdivErrorEstimate.cpp TPZMixedHdivErrorEstimate - TPZMatLaplacianHybrid.cpp - TPZMatLaplacianHybrid.h - TPZPressureProjection.h - TPZPressureProjection.cpp - TPZHybridH1ErrorEstimateMaterial.h - TPZHybridH1ErrorEstimateMaterial.cpp - TPZSteklovMaterial.h - TPZSteklovMaterial.cpp + #TPZMatLaplacianHybrid.cpp + #TPZMatLaplacianHybrid.h + #TPZPressureProjection.h + #TPZPressureProjection.cpp + #TPZHybridH1ErrorEstimateMaterial.h + #TPZHybridH1ErrorEstimateMaterial.cpp + #TPZSteklovMaterial.h + #TPZSteklovMaterial.cpp ) target_include_directories(ErrorEstimationLib PUBLIC ${CMAKE_CURRENT_LIST_DIR}) diff --git a/ErrorEstimation/Material/TPZHDivErrorEstimateMaterial.cpp b/ErrorEstimation/Material/TPZHDivErrorEstimateMaterial.cpp index 6147cc99..e827a48f 100644 --- a/ErrorEstimation/Material/TPZHDivErrorEstimateMaterial.cpp +++ b/ErrorEstimation/Material/TPZHDivErrorEstimateMaterial.cpp @@ -8,53 +8,37 @@ #include "TPZHDivErrorEstimateMaterial.h" #include "pzaxestools.h" #include "TPZAnalyticSolution.h" -#include "TPZMaterialDataT.h" - - +#include "TPZBndCond.h" #ifdef LOG4CXX static LoggerPtr logger(Logger::getLogger("pz.errorestimation.hdiv")); #endif -TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(int matid, int dim) : TPZMixedPoisson(matid,dim) -{ - -} +TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(int matid, int dim) : TPZMixedDarcyFlow(matid, dim) {} -TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial() : TPZMixedPoisson() -{ - -} +TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial() : TPZMixedDarcyFlow() {} -TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(const TPZMixedPoisson ©) : TPZMixedPoisson(copy) -{ - -} +TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(const TPZMixedDarcyFlow ©) + : TPZMixedDarcyFlow(copy) {} -TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(const TPZHDivErrorEstimateMaterial ©) : TPZMixedPoisson(copy) -{ - -} +TPZHDivErrorEstimateMaterial::TPZHDivErrorEstimateMaterial(const TPZHDivErrorEstimateMaterial ©) + : TPZMixedDarcyFlow(copy) {} -TPZHDivErrorEstimateMaterial::~TPZHDivErrorEstimateMaterial() -{ - -} +TPZHDivErrorEstimateMaterial::~TPZHDivErrorEstimateMaterial() = default; -TPZHDivErrorEstimateMaterial &TPZHDivErrorEstimateMaterial::operator=(const TPZHDivErrorEstimateMaterial ©) -{ - TPZMixedPoisson::operator=(copy); +TPZHDivErrorEstimateMaterial &TPZHDivErrorEstimateMaterial::operator=(const TPZHDivErrorEstimateMaterial ©) { + TPZMixedDarcyFlow::operator=(copy); return *this; } -int TPZHDivErrorEstimateMaterial::FirstNonNullApproxSpaceIndex(const TPZVec> &datavec) const { +int TPZHDivErrorEstimateMaterial::FirstNonNullApproxSpaceIndex(const TPZVec> &datavec) { int nvec = datavec.NElements(); int firstNoNullposition = -1; - for(int ivec = 0; ivec < nvec ; ivec++){ - auto shapetype = datavec[ivec].fShapeType; - if(shapetype != TPZMaterialData::EEmpty){ + for (int ivec = 0; ivec < nvec; ivec++) { + TPZMaterialData::MShapeFunctionType shapetype = datavec[ivec].fShapeType; + if (shapetype != TPZMaterialData::EEmpty) { firstNoNullposition = ivec; return firstNoNullposition; } @@ -63,8 +47,8 @@ int TPZHDivErrorEstimateMaterial::FirstNonNullApproxSpaceIndex(const TPZVec> &datavec, REAL weight, TPZFMatrix &ek, TPZFMatrix &ef) -{ +void TPZHDivErrorEstimateMaterial::Contribute(const TPZVec> &datavec, REAL weight, + TPZFMatrix &ek, TPZFMatrix &ef) { /** datavec[0] H1 mesh, local uk/grad v for Mark reconstruction and Empty for H1 reconstruction @@ -81,114 +65,94 @@ void TPZHDivErrorEstimateMaterial::Contribute(const TPZVec &phiuk = datavec[H1functionposition].phi; TPZFMatrix &dphiukaxes = datavec[H1functionposition].dphix; - TPZFNMatrix<9,REAL> dphiuk(3,dphiukaxes.Cols()); + TPZFNMatrix<9, REAL> dphiuk(3, 1); TPZAxesTools::Axes2XYZ(dphiukaxes, dphiuk, datavec[H1functionposition].axes); - - - int nphiuk = phiuk.Rows(); - - TPZFMatrix solsigmafem(3,nphiuk),solukfem(1,1); + + TPZFMatrix solsigmafem(3, 1), solukfem(1, 1); solsigmafem.Zero(); solukfem.Zero(); - - - - //potetial fem - solukfem(0,0) = datavec[3].sol[0][0]; - //flux fem - for (int ip = 0; ip<3; ip++){ + //potetial fem + solukfem(0, 0) = datavec[3].sol[0][0]; + //flux fem + for (int ip = 0; ip < 3; ip++) { - solsigmafem(ip,0) = datavec[2].sol[0][ip]; - } + solsigmafem(ip, 0) = datavec[2].sol[0][ip]; + } - - - TPZFNMatrix<9,REAL> PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); STATE perm = GetPermeability(datavec[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); - - - TPZFMatrix kgraduk(3,nphiuk,0.); - - - for(int irow=0 ; irow kgraduk(3, nphiuk, 0.); + + for (int irow = 0; irow < nphiuk; irow++) { + //K graduk - for(int id=0; id< dim; id++){ - - for(int jd=0; jd< dim;jd++){ - - kgraduk(id,irow) += PermTensor(id,jd)*dphiuk(jd,irow); - - } + for (int id = 0; id < dim; id++) { + + kgraduk(id, irow) += perm * dphiuk(id, irow); //bk = (-1)*int_k sigmaukfem.grad phi_i,here dphiuk is multiplied by axes //the minus sign is necessary because we are working with sigma_h = - K grad u, Mark works with sigma_h = K grad u - - ef(irow,0)+=(-1.)*weight*dphiuk(id,irow)*solsigmafem(id,0); + + ef(irow, 0) += (-1.) * weight * dphiuk(id, irow) * solsigmafem(id, 0); } - + //matrix Sk= int_{K} K graduk.gradv - for(int jcol=0; jcol> &datavec, REAL weight, TPZFMatrix &ek, - TPZFMatrix &ef, TPZBndCondT &bc) { +void TPZHDivErrorEstimateMaterial::ContributeBC(const TPZVec> &datavec, REAL weight, + TPZFMatrix &ek, TPZFMatrix &ef, TPZBndCondT &bc) { /* Add Robin boundary condition for local problem ek+= ef+= */ - int H1functionposition = FirstNonNullApproxSpaceIndex(datavec); - int dim = datavec[H1functionposition].axes.Rows(); + auto H1functionposition = FirstNonNullApproxSpaceIndex(datavec); + auto dim = datavec[H1functionposition].axes.Rows(); TPZFMatrix &phi_i = datavec[H1functionposition].phi; - int nphi_i = phi_i.Rows(); + auto nphi_i = phi_i.Rows(); TPZFMatrix solsigmafem(3, 1); solsigmafem.Zero(); @@ -200,11 +164,6 @@ void TPZHDivErrorEstimateMaterial::ContributeBC(const TPZVec PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); - STATE perm = GetPermeability(datavec[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); if (bc.HasForcingFunctionBC()) { TPZManVector res(3); @@ -212,10 +171,11 @@ void TPZHDivErrorEstimateMaterial::ContributeBC(const TPZVec > &datavec) const { +void TPZHDivErrorEstimateMaterial::FillDataRequirements(TPZVec> &datavec) const { - //fem solution for flux and potential - datavec[2].SetAllRequirements(false); - datavec[2].fNeedsSol = true; - datavec[2].fNeedsNormal = true; - - datavec[3].SetAllRequirements(false); - datavec[3].fNeedsSol = true; - + //fem solution for flux and potential + datavec[2].SetAllRequirements(false); + datavec[2].fNeedsSol = true; + datavec[2].fNeedsNormal = true; + datavec[3].SetAllRequirements(false); + datavec[3].fNeedsSol = true; } -void TPZHDivErrorEstimateMaterial::FillBoundaryConditionDataRequirements(int type,TPZVec > &datavec) const { - +void +TPZHDivErrorEstimateMaterial::FillBoundaryConditionDataRequirements(int type, + TPZVec> &datavec) const { datavec[2].SetAllRequirements(false); datavec[2].fNeedsSol = true; datavec[2].fNeedsNormal = true; - - } -void TPZHDivErrorEstimateMaterial::Errors(const TPZVec> &data, TPZVec &errors) -{ - /** - datavec[0] H1 mesh, uh_reconstructed - datavec[1] L2 mesh, - datavec[2] Hdiv fem mesh, sigma_h - datavec[3] L2 mesh fem, u_h - - error[0] - error computed with exact pressure - error[1] - error computed with reconstructed pressure - error[2] - energy error computed with exact solution - error[3] - energy error computed with reconstructed solution - error[4] - oscilatory data error - **/ - - errors.Resize(NEvalErrors()); - errors.Fill(0.0); - - - TPZManVector fluxfem(3); - STATE divsigmafem, pressurefem, pressurereconstructed; - divsigmafem = 0.; - pressurefem = 0.; - pressurereconstructed = 0.; - - - TPZFNMatrix<3,REAL> fluxreconstructed(3,1), fluxreconstructed2(3,1); - - - fluxfem = data[2].sol[0]; - - - divsigmafem= data[2].divsol[0][0]; - - STATE divtest=0.; - for (int j=0; j> &data, TPZVec &errors) { + /* + * data[0]: H1 mesh, uh_reconstructed + * data[1]: L2 mesh, + * data[2]: Hdiv fem mesh, sigma_h + * data[3]: L2 mesh fem, u_h + * + * error[0]: error computed with exact/reference pressure + * error[1]: error computed with reconstructed pressure + * error[2]: energy error computed with exact/reference solution + * error[3]: energy error computed with reconstructed solution + * error[4]: oscillatory data error + */ + + // Variables to store error norms/components + STATE flux_error_estimate = 0.; + STATE flux_error_exact = 0.; + STATE pressure_error_estimate = 0.; + STATE pressure_error_exact = 0.; + STATE residual_error = 0.; + + const auto perm = GetPermeability(data[1].x); + + const auto pressurefem = data[3].sol[0][0]; + const auto &fluxfem = data[2].sol[0]; + const auto divfluxfem = data[2].divsol[0][0]; + const auto pressurerec = data[1].sol[0][0]; + + // Calculate estimated errors + const auto &du_rec = data[1].dsol[0]; + TPZFNMatrix<3, REAL> flux_rec(3, 1); + TPZAxesTools::Axes2XYZ(du_rec, flux_rec, data[1].axes); + pressure_error_estimate = (pressurefem - pressurerec) * (pressurefem - pressurerec); + for (int i = 0; i < 3; i++) { + flux_rec(i, 0) = -perm * flux_rec(i, 0); + flux_error_estimate += (fluxfem[i] - flux_rec(i, 0)) * (fluxfem[i] - flux_rec(i, 0)) / perm; } - int H1functionposition = 0; - H1functionposition = FirstNonNullApproxSpaceIndex(data); - - - TPZVec divsigma(1); - divsigma[0]=0.; - - TPZManVector u_exact(1); - TPZFNMatrix<9, STATE> du_exact(3, 3); - if(this->fExactSol){ - - this->fExactSol(data[H1functionposition].x,u_exact,du_exact); - - } - if(this->HasForcingFunction()) - { - this->ForcingFunction()(data[H1functionposition].x,divsigma); + // Calculate exact errors, if applicable + if (this->fExactSol) { + TPZManVector u_exact(1); + TPZFNMatrix<3, STATE> flux_exact(3, 1); + this->fExactSol(data[1].x, u_exact, flux_exact); + pressure_error_exact = (pressurefem - u_exact[0]) * (pressurefem - u_exact[0]); + for (int i = 0; i < 3; i++) { + flux_exact(i, 0) = -perm * flux_exact(i, 0); + flux_error_exact += (fluxfem[i] - flux_exact(i, 0)) * (fluxfem[i] - flux_exact(i, 0)) / perm; + } + } else if (this->fHasReferenceSolution) { - } - - REAL residual = 0.; + const auto u_ref = data[5].sol[0][0]; + const auto &flux_ref = data[4].sol[0]; - residual = (divsigma[0] - divsigmafem)*(divsigma[0] - divsigmafem); - - - pressurereconstructed = data[H1functionposition].sol[0][0]; - - - pressurefem = data[3].sol[0][0]; - - - TPZFNMatrix<9,REAL> PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); - STATE perm = GetPermeability(data[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); - - TPZFNMatrix<3,REAL> fluxexactneg; - - //sigmarec = -K grad(urec) - // sigmak = -K graduk - - { - TPZFNMatrix<9,REAL> gradpressure(3,1); - for (int i=0; i<3; i++) { - gradpressure(i,0) = du_exact[i]; - } - PermTensor.Multiply(gradpressure,fluxexactneg); - } - - - - TPZFMatrix &dsolaxes = data[H1functionposition].dsol[0]; - TPZFNMatrix<9,REAL> fluxrec(3,0); - TPZAxesTools::Axes2XYZ(dsolaxes, fluxrec, data[H1functionposition].axes); - - for(int id=0 ; id<3; id++) { - fluxreconstructed2(id,0) = (-1.)*fluxrec(id,0); - } - - PermTensor.Multiply(fluxreconstructed2,fluxreconstructed); - - - REAL innerexact = 0.; - REAL innerestimate = 0.; - - - -#ifdef ERRORESTIMATION_DEBUG2 - std::cout<<"flux fem "<isDebugEnabled()) { - std::stringstream sout; - sout << "Coord: " << data[H1functionposition].x[0] << ", " << data[H1functionposition].x[1] << ", " - << data[H1functionposition].x[2] << '\n'; - sout << "PressureReconstructed = " << pressurereconstructed << "\n"; - sout << "FluxReconstructed = " << fluxreconstructed[0] << ", " << fluxreconstructed[1] << ", " - << fluxreconstructed[2] << "\n"; - LOGPZ_DEBUG(logger, sout.str()) + // Calculate residual component + TPZManVector source_term(1, 0); + if (this->fForcingFunction) { + this->fForcingFunction(data[1].x, source_term); } -#endif + residual_error = (source_term[0] - divfluxfem) * (source_term[0] - divfluxfem); + + // Fill error vector + errors[0] = pressure_error_exact; + errors[1] = pressure_error_estimate; + errors[2] = flux_error_exact; + errors[3] = flux_error_estimate; + errors[4] = residual_error; } -int TPZHDivErrorEstimateMaterial::VariableIndex(const std::string &name) const -{ - if(name == "FluxFem") return 40; - if(name == "FluxReconstructed") return 41; - if(name == "FluxExact") return 42; - if(name == "PressureFem") return 43; - if(name == "PressureReconstructed") return 44; - if(name == "PressureExact") return 45; - if(name == "PressureErrorExact") return 100; - if(name == "PressureErrorEstimate") return 101; - if(name == "EnergyErrorExact") return 102; - if(name == "EnergyErrorEstimate") return 103; - if(name == "ResidualError") return 104; - if(name == "PressureEffectivityIndex") return 105; - if(name == "EnergyEffectivityIndex") return 106; - if(name == "POrder") return 46; - +int TPZHDivErrorEstimateMaterial::VariableIndex(const std::string &name) const { + if (name == "FluxFem") return 40; + if (name == "FluxReconstructed") return 41; + if (name == "FluxExact") return 42; + if (name == "PressureFem") return 43; + if (name == "PressureReconstructed") return 44; + if (name == "PressureExact") return 45; + if (name == "POrder") return 46; + if (name == "Permeability") return 47; + if (name == "PressureErrorExact") return 100; + if (name == "PressureErrorEstimate") return 101; + if (name == "EnergyErrorExact") return 102; + if (name == "EnergyErrorEstimate") return 103; + if (name == "ResidualError") return 104; + if (name == "PressureEffectivityIndex") return 105; + if (name == "EnergyEffectivityIndex") return 106; + return -1; } -int TPZHDivErrorEstimateMaterial::NSolutionVariables(int var) const -{ +int TPZHDivErrorEstimateMaterial::NSolutionVariables(int var) const { switch (var) { case 40: case 41: case 42: return 3; - break; case 43: case 44: case 45: @@ -467,7 +352,6 @@ int TPZHDivErrorEstimateMaterial::NSolutionVariables(int var) const case 105: case 106: return 1; - break; default: DebugStop(); break; @@ -482,8 +366,8 @@ int TPZHDivErrorEstimateMaterial::NSolutionVariables(int var) const * @param Solout [out] is the solution vector */ -void TPZHDivErrorEstimateMaterial::Solution(const TPZVec> &datavec, int var, TPZVec &Solout) -{ +void +TPZHDivErrorEstimateMaterial::Solution(const TPZVec> &datavec, int var, TPZVec &Solout) { /** datavec[0] H1 mesh, uh_reconstructed for Mark reconstruction and Empty for H1 reconstruction @@ -491,58 +375,47 @@ void TPZHDivErrorEstimateMaterial::Solution(const TPZVec datavec[2] Hdiv fem mesh, sigma_h datavec[3] L2 mesh fem, u_h **/ - - int H1functionposition = 0; - H1functionposition = FirstNonNullApproxSpaceIndex(datavec); - - TPZFNMatrix<9,REAL> PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); + + int H1functionposition = FirstNonNullApproxSpaceIndex(datavec); + STATE perm = GetPermeability(datavec[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); - - - - TPZManVector pressexact(1,0.); - TPZFNMatrix<9,STATE> gradu(3,1,0.), fluxinv(3,1); - - if(HasExactSol()) - { - this->ExactSol()(datavec[H1functionposition].x, pressexact,gradu); - + + TPZManVector pressexact(1, 0.); + TPZFNMatrix<9, STATE> gradu(3, 1, 0.), fluxinv(3, 1); + + if (fExactSol) { + this->fExactSol(datavec[H1functionposition].x, pressexact, gradu); + } + + for (int i = 0; i < 3; i++) { + fluxinv(i, 0) = perm * gradu(i, 0); } - - PermTensor.Multiply(gradu, fluxinv); - - int dim=this->fDim; - switch (var) - { + + int dim = this->fDim; + switch (var) { case 40://FluxFem - - for(int i=0; i<3; i++) Solout[i] = datavec[2].sol[0][i]; + + for (int i = 0; i < 3; i++) Solout[i] = datavec[2].sol[0][i]; break; - case 41:{//FluxReconstructed is grad U + case 41: {//FluxReconstructed is grad U TPZFMatrix &dsolaxes = datavec[H1functionposition].dsol[0]; - TPZFNMatrix<9,REAL> dsol(3,0); - TPZFNMatrix<9,REAL> KGradsol(3,0); + TPZFNMatrix<9, REAL> dsol(3, 1); TPZAxesTools::Axes2XYZ(dsolaxes, dsol, datavec[H1functionposition].axes); - - PermTensor.Multiply(dsol,KGradsol); - - for(int id=0 ; id case 46://order p Solout[0] = datavec[1].p; break; + case 47: // Permeability + Solout[0] = perm; + break; default: DebugStop(); } } - -void TPZHDivErrorEstimateMaterial:: ErrorsBC(TPZVec> &data, TPZVec &errors, TPZBndCondT &bc){ - - if(bc.Type()== 4){ - int H1functionposition = 0; - H1functionposition = FirstNonNullApproxSpaceIndex(data); - TPZVec u_exact(1); - TPZFMatrix du_exact(3,1,0.); - if(HasExactSol()) - { - ExactSol()(data[H1functionposition].x,u_exact,du_exact); - } - else - { - return; - } - - - errors.Resize(NEvalErrors()); - errors.Fill(0.0); - - - TPZFNMatrix<3,REAL> fluxreconstructed(3,1), fluxreconstructed2(3,1); - TPZManVector fluxfem(3); - - - REAL normalsigmafem = 0.,normalsigmarec = 0.,urec=0.;; - normalsigmafem = data[2].sol[0][0];// sigma.n - urec = data[H1functionposition].sol[0][0]; - - - - REAL u_D = 0.,g = 0.; - REAL normflux = 0.; - - TPZManVector fluxrec(fDim); - this->Solution(data,VariableIndex("FluxReconstructed"), fluxrec); - - // std::cout<<"flux_rec "< PermTensor, InvPermTensor; - TPZManVector res(3); - TPZFNMatrix<9, STATE> gradu(this->Dimension(), 1); - u_D = u_exact[0]; - std::cout << "MUDEI O FLUXO DO PROGRAMA OLHE COM CUIDADO\n"; - - - for(int i=0; i<3; i++) - { - for(int j=0; j<3; j++) - { - - normflux += data[2].normal[i]*PermTensor(i,j)*gradu(j,0); - - } - } - g = (-1)*normflux; - - // std::cout<<"n_0 "< +#include #include "DarcyFlow/TPZMixedDarcyFlow.h" +#include "TPZMaterialDataT.h" -typedef TPZMixedDarcyFlow TPZMixedPoisson; - -class TPZHDivErrorEstimateMaterial : public TPZMixedPoisson { +class TPZHDivErrorEstimateMaterial: public TPZMixedDarcyFlow { public: TPZHDivErrorEstimateMaterial(int matid, int dim); @@ -25,20 +24,23 @@ class TPZHDivErrorEstimateMaterial : public TPZMixedPoisson { TPZHDivErrorEstimateMaterial(const TPZHDivErrorEstimateMaterial ©); - TPZHDivErrorEstimateMaterial(const TPZMixedPoisson ©); + explicit TPZHDivErrorEstimateMaterial(const TPZMixedDarcyFlow ©); - virtual ~TPZHDivErrorEstimateMaterial(); + ~TPZHDivErrorEstimateMaterial() override; TPZHDivErrorEstimateMaterial &operator=(const TPZHDivErrorEstimateMaterial ©); - virtual void Contribute(const TPZVec> &datavec, REAL weight, TPZFMatrix &ek, + void Contribute(const TPZVec> &datavec, REAL weight, TPZFMatrix &ek, TPZFMatrix &ef) override; - virtual void ContributeBC(const TPZVec> &datavec, REAL weight, TPZFMatrix &ek, - TPZFMatrix &ef, TPZBndCondT &bc) override; + void ContributeBC(const TPZVec> &datavec, REAL weight, TPZFMatrix &ek, + TPZFMatrix &ef, TPZBndCondT &bc) override; + + // TODO add doc + void FillDataRequirements(TPZVec > &datavec) const override; - virtual void FillDataRequirements(TPZVec> &datavec) const override; - virtual void FillBoundaryConditionDataRequirements(int type, TPZVec> &datavec) const override; + // TODO add doc + void FillBoundaryConditionDataRequirements(int type, TPZVec > &datavec) const override; bool fNeumannLocalProblem = false; @@ -46,25 +48,33 @@ class TPZHDivErrorEstimateMaterial : public TPZMixedPoisson { fNeumannLocalProblem = neumannProblem; } - virtual int NEvalErrors() const override { return 5; } // erro de oscilacao de dados tbem + [[nodiscard]] int NEvalErrors() const override { return 5; } /// Compute the error and error estimate // error[0] - error computed with exact pressure // error[1] - error computed with reconstructed pressure // error[2] - energy error computed with exact solution // error[3] - energy error computed with reconstructed solution - virtual void Errors(const TPZVec> &data, TPZVec &errors) override; + void Errors(const TPZVec> &data, TPZVec &errors) override; + + //void ErrorsBC(const TPZVec> &data, TPZVec &u_exact, TPZFMatrix &du_exact, + // TPZVec &errors, TPZBndCond &bc); - void ErrorsBC(TPZVec> &data, - TPZVec &errors, TPZBndCondT &bc); + [[nodiscard]] int VariableIndex(const std::string &name) const override; + [[nodiscard]] int NSolutionVariables(int var) const override; - virtual int VariableIndex(const std::string &name) const override; - virtual int NSolutionVariables(int var) const override; - virtual void Solution(const TPZVec> &datavec, int var, TPZVec &Solout) override; + void Solution(const TPZVec> &datavec, int var, TPZVec &Solout) override; // Returns the first non-null approximation space index, which will be the // H1 reconstruction space - int FirstNonNullApproxSpaceIndex(const TPZVec> &datavec) const; + static int FirstNonNullApproxSpaceIndex(const TPZVec> &datavec); + + void SetReferenceSolution(bool hasRefSol) { + fHasReferenceSolution = hasRefSol; + } + +private: + bool fHasReferenceSolution{false}; }; #endif /* TPZHDivErrorEstimateMaterial_hpp */ diff --git a/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.cpp b/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.cpp index 4d1b4a55..e802b024 100644 --- a/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.cpp +++ b/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.cpp @@ -6,64 +6,44 @@ // #include "TPZMixedHdivErrorEstimate.h" -#include "DarcyFlow/TPZMixedDarcyFlow.h" #include "pzaxestools.h" -#include "TPZMaterialDataT.h" - -template -TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate() : MixedMat() +TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate() : TPZMixedDarcyFlow() { } -template -TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(int matid, int dim) : MixedMat(matid,dim) +TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(int matid, int dim) : TPZMixedDarcyFlow(matid,dim) { } -template -TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(const MixedMat ©) : MixedMat(copy) +TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(const TPZMixedDarcyFlow ©) : TPZMixedDarcyFlow(copy) { } -template -TPZMixedHDivErrorEstimate::~TPZMixedHDivErrorEstimate() -{ - -} +TPZMixedHDivErrorEstimate::~TPZMixedHDivErrorEstimate() = default; -template -TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(const TPZMixedHDivErrorEstimate &cp) : MixedMat(cp) +TPZMixedHDivErrorEstimate::TPZMixedHDivErrorEstimate(const TPZMixedHDivErrorEstimate &cp) : TPZMixedDarcyFlow(cp) { } -template -TPZMixedHDivErrorEstimate &TPZMixedHDivErrorEstimate::operator=(const TPZMixedHDivErrorEstimate ©) +TPZMixedHDivErrorEstimate &TPZMixedHDivErrorEstimate::operator=(const TPZMixedHDivErrorEstimate ©) { - MixedMat *me = this; - me->operator=(copy); + TPZMixedDarcyFlow::operator=(copy); return *this; } -template -void TPZMixedHDivErrorEstimate::FillDataRequirements(TPZVec > &datavec) const -{ - MixedMat::FillDataRequirements(datavec); - { - int i = 1; - datavec[i].SetAllRequirements(false); - //datavec[0].fNeedsSol = true; - datavec[i].fNeedsSol = true; - } - +void TPZMixedHDivErrorEstimate::FillDataRequirements(TPZVec> &datavec) const { + TPZMixedDarcyFlow::FillDataRequirements(datavec); + int i = 1; + datavec[i].SetAllRequirements(false); + datavec[i].fNeedsSol = true; } -template -int TPZMixedHDivErrorEstimate::VariableIndex(const std::string &name) const +int TPZMixedHDivErrorEstimate::VariableIndex(const std::string &name) const { if(name == "FluxFem") return 40; if(name == "FluxReconstructed") return 41; @@ -83,15 +63,12 @@ int TPZMixedHDivErrorEstimate::VariableIndex(const std::string &name) return -1; } -template -int TPZMixedHDivErrorEstimate::NSolutionVariables(int var) const -{ +int TPZMixedHDivErrorEstimate::NSolutionVariables(int var) const { switch (var) { case 40: case 41: case 42: return 3; - break; case 43: case 44: case 45: @@ -104,7 +81,6 @@ int TPZMixedHDivErrorEstimate::NSolutionVariables(int var) const case 105: case 106: return 1; - break; default: DebugStop(); break; @@ -118,9 +94,8 @@ int TPZMixedHDivErrorEstimate::NSolutionVariables(int var) const * @param var [in] number of solution variables. See NSolutionVariables() method * @param Solout [out] is the solution vector */ -template -void TPZMixedHDivErrorEstimate::Solution(const TPZVec> &datavec, int var, TPZVec &Solout) -{ +void +TPZMixedHDivErrorEstimate::Solution(const TPZVec> &datavec, int var, TPZVec &Solout) { /** datavec[0]= Hdiv Resconstructed datavec[1]= Pressure Resconstructed @@ -128,30 +103,22 @@ void TPZMixedHDivErrorEstimate::Solution(const TPZVec PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); - - STATE perm = MixedMat::GetPermeability(datavec[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); - int dim = MixedMat::fDim; + const auto perm = GetPermeability(datavec[1].x); - STATE pressureexact = 0.; - TPZManVector pressvec(1,0.); - TPZFNMatrix<9, STATE> gradu(dim, 1, 0.), fluxinv(dim, 1); + TPZManVector pressexact(1, 0.); + TPZFNMatrix<9, STATE> gradu(3, 1, 0.), fluxinv(3, 1); - if(MixedMat::fExactSol) - { + if (fExactSol) { + this->fExactSol(datavec[0].x, pressexact, gradu); + } - MixedMat::fExactSol(datavec[0].x, pressvec,gradu); - gradu.Resize(3, 1); - //gradu(2,0) = 0.; + for (int i = 0; i < 3; i++) { + fluxinv(i, 0) = perm * gradu(i, 0); } - - PermTensor.Multiply(gradu, fluxinv); - pressureexact = pressvec[0]; + + int dim = this->fDim; + switch (var) { case 40://FluxFem @@ -170,7 +137,7 @@ void TPZMixedHDivErrorEstimate::Solution(const TPZVec::Solution(const TPZVec -void TPZMixedHDivErrorEstimate::Errors(const TPZVec > &data, TPZVec &errors) -{ - +void TPZMixedHDivErrorEstimate::Errors(const TPZVec> &data, TPZVec &errors) { + /** datavec[0]= Hdiv Reconstructed datavec[1]= Pressure Reconstructed datavec[2]= Hdiv FEM datavec[3]= Pressure FEM **/ - + errors.Resize(NEvalErrors()); errors.Fill(0.0); - const int dim = MixedMat::fDim; + const int dim = this->fDim; TPZManVector fluxfem(3), fluxreconstructed(3), pressurefem(1), pressurereconstructed(1); - + fluxreconstructed = data[0].sol[0]; fluxfem = data[2].sol[0]; STATE divsigmafem=data[2].divsol[0][0]; - - + + TPZVec divsigma(1,0.); TPZManVector u_exact(1); TPZFNMatrix<9, STATE> du_exact(3, 3); - if(this->HasExactSol()){ - this->ExactSol()(data[0].x,u_exact,du_exact); + if(this->fExactSol){ + this->fExactSol(data[0].x,u_exact,du_exact); + this->fForcingFunction(data[0].x,divsigma); } - if(this->ForcingFunction()) - { - this->ForcingFunction()(data[0].x,divsigma); - } - - REAL residual = 0.; -// std::cout<<" divsigma "< PermTensor(3,3); - TPZFNMatrix<9,REAL> InvPermTensor(3,3); - - REAL perm = MixedMat::GetPermeability(data[1].x); - PermTensor.Diagonal(perm); - InvPermTensor.Diagonal(1./perm); - - - TPZFNMatrix<3,REAL> fluxexactneg; + + const auto perm = GetPermeability(data[1].x); //sigma=-K grad(u) + TPZFNMatrix<3, REAL> fluxexactneg(3, 1); + { - TPZFNMatrix<9,REAL> gradpressure(dim,1); - for (int i=0; i gradpressure(3, 1); + for (int i = 0; i < 3; i++) { + fluxexactneg(i, 0) = perm * du_exact[i]; } - PermTensor.Multiply(gradpressure,fluxexactneg); } - - + REAL innerexact = 0.; REAL innerestimate = 0.; - for (int i=0; i::Errors(const TPZVec; diff --git a/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.h b/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.h index b920d295..1186b939 100644 --- a/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.h +++ b/ErrorEstimation/Material/TPZMixedHdivErrorEstimate.h @@ -8,58 +8,43 @@ #ifndef TPZMixedErrorEstimate_hpp #define TPZMixedErrorEstimate_hpp -#include +#include +#include +#include #include "pzreal.h" -template -class TPZFMatrix; -class TPZMaterial; -class TPZMaterialData; +class TPZMixedHDivErrorEstimate : public TPZMixedDarcyFlow { -template -class TPZMaterialDataT; - -template -class TPZVec; - -template -class TPZMixedHDivErrorEstimate : public MixedMat -{ - - public: TPZMixedHDivErrorEstimate(); TPZMixedHDivErrorEstimate(int matid, int dim); - virtual ~TPZMixedHDivErrorEstimate(); + ~TPZMixedHDivErrorEstimate() override; - TPZMixedHDivErrorEstimate(const MixedMat &cp); + explicit TPZMixedHDivErrorEstimate(const TPZMixedDarcyFlow &cp); - TPZMixedHDivErrorEstimate(const TPZMixedHDivErrorEstimate &cp); + TPZMixedHDivErrorEstimate(const TPZMixedHDivErrorEstimate &cp); TPZMixedHDivErrorEstimate &operator=(const TPZMixedHDivErrorEstimate ©); - virtual TPZMaterial * NewMaterial() const override { + [[nodiscard]] TPZMaterial * NewMaterial() const override { return new TPZMixedHDivErrorEstimate(*this); } - void FillDataRequirements(TPZVec > &datavec) const override; + void FillDataRequirements(TPZVec> &datavec) const override; /// make a contribution to the error computation - virtual void Errors(const TPZVec> &data, TPZVec &errors) override; - - virtual int NEvalErrors() const override { - return 5; - - } + void Errors(const TPZVec> &data, TPZVec &errors) override; + + [[nodiscard]] int NEvalErrors() const override { return 5; } - virtual int VariableIndex(const std::string &name) const override; + [[nodiscard]] int VariableIndex(const std::string &name) const override; - virtual int NSolutionVariables(int var) const override; + [[nodiscard]] int NSolutionVariables(int var) const override; /** * @brief It return a solution to multiphysics simulation. @@ -67,7 +52,7 @@ class TPZMixedHDivErrorEstimate : public MixedMat * @param var [in] number of solution variables. See NSolutionVariables() method * @param Solout [out] is the solution vector */ - virtual void Solution(const TPZVec> &datavec, int var, TPZVec &Solout) override; + void Solution(const TPZVec> &datavec, int var, TPZVec &Solout) override; /// whether the post processing mesh will be H(div) or H1 bool fPostProcesswithHDiv = true; diff --git a/ErrorEstimation/ProblemConfig.h b/ErrorEstimation/ProblemConfig.h index 229d69b0..1b53cb47 100644 --- a/ErrorEstimation/ProblemConfig.h +++ b/ErrorEstimation/ProblemConfig.h @@ -14,9 +14,9 @@ /// class to guide the error estimator struct ProblemConfig { - + TPZCompMesh *cmesh = nullptr; /// geometric mesh on which the computational meshes are based - TPZGeoMesh *gmesh = 0; + TPZGeoMesh *gmesh = nullptr; /// polynomial order of the original mesh int porder = 2; /// increment in internal order of flux and pressure @@ -41,7 +41,7 @@ struct ProblemConfig bool GalvisExample = false; bool TensorNonConst = false; bool MeshNonConvex = false; - + STATE alpha=1; STATE Km = 0.; STATE coefG = 0.; diff --git a/ErrorEstimation/TPZCreateMultiphysicsSpace.cpp b/ErrorEstimation/TPZCreateMultiphysicsSpace.cpp index 62a86ae1..ab9c5af5 100644 --- a/ErrorEstimation/TPZCreateMultiphysicsSpace.cpp +++ b/ErrorEstimation/TPZCreateMultiphysicsSpace.cpp @@ -14,7 +14,7 @@ #include "pzelementgroup.h" #include "pzcondensedcompel.h" #include "TPZNullMaterial.h" -#include "TPZLagrangeMultiplierCS.h" +#include "TPZLagrangeMultiplier.h" #include "TPZMultiphysicsCompMesh.h" #include "TPZMultiphysicsInterfaceEl.h" #include "TPZHybridizeHDiv.h" @@ -270,7 +270,7 @@ void TPZCreateMultiphysicsSpace::CreatePressureBoundaryElements(TPZCompMesh *pre } pressure->ExpandSolution(); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { int err = 0; std::map gelindexes; @@ -312,7 +312,7 @@ void TPZCreateMultiphysicsSpace::CreatePressureBoundaryElements(TPZCompMesh *pre void TPZCreateMultiphysicsSpace::InsertPressureMaterialIds(TPZCompMesh *pressure) { for (auto matid:fMaterialIds) { - TPZNullMaterial *nullmat = new TPZNullMaterial(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); pressure->InsertMaterialObject(nullmat); @@ -320,21 +320,21 @@ void TPZCreateMultiphysicsSpace::InsertPressureMaterialIds(TPZCompMesh *pressure if((fH1Hybrid.fHybridizeBCLevel == 2) || (fH1Hybrid.fHybridizeBCLevel == 0)) { for (auto matid:fBCMaterialIds) { - TPZNullMaterial *nullmat = new TPZNullMaterial<>(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); pressure->InsertMaterialObject(nullmat); } } { - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(fH1Hybrid.fMatWrapId); + auto *nullmat = new TPZNullMaterial(fH1Hybrid.fMatWrapId); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); pressure->InsertMaterialObject(nullmat); } if(fSpaceType == EH1HybridSquared) { - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(fH1Hybrid.fInterfacePressure); + auto *nullmat = new TPZNullMaterial(fH1Hybrid.fInterfacePressure); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); pressure->InsertMaterialObject(nullmat); @@ -347,7 +347,7 @@ void TPZCreateMultiphysicsSpace::InsertFluxMaterialIds(TPZCompMesh *fluxmesh) { if (fSpaceType == EH1Hybrid || fSpaceType == EH1HybridSquared) { int matid = fH1Hybrid.fFluxMatId; - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(fDimension-1); nullmat->SetNStateVariables(1); fluxmesh->InsertMaterialObject(nullmat); @@ -358,7 +358,7 @@ void TPZCreateMultiphysicsSpace::InsertFluxMaterialIds(TPZCompMesh *fluxmesh) if(fH1Hybrid.fHybridizeBCLevel == 1) { for (auto matid:fBCMaterialIds) { - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); fluxmesh->InsertMaterialObject(nullmat); @@ -370,7 +370,7 @@ void TPZCreateMultiphysicsSpace::InsertFluxMaterialIds(TPZCompMesh *fluxmesh) void TPZCreateMultiphysicsSpace::InsertNullSpaceMaterialIds(TPZCompMesh *nullspace) { for (auto matid:fMaterialIds) { - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(fDimension); nullmat->SetNStateVariables(1); nullspace->InsertMaterialObject(nullmat); @@ -524,7 +524,7 @@ void TPZCreateMultiphysicsSpace::AddInterfaceElements(TPZMultiphysicsCompMesh *m fluxcandidate = gelsideflux.HasLowerLevelNeighbour(fH1Hybrid.fFluxMatId); DebugStop(); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG if(fluxcandidate == gelsideflux) { DebugStop(); @@ -659,7 +659,7 @@ TPZCompEl *TPZCreateMultiphysicsSpace::FindFluxElement(TPZCompEl *wrapelement) } } } - return NULL; + return nullptr; } /// Compute Periferal Material ids @@ -691,7 +691,7 @@ static void InsertNullMaterial(int matid, int dim, int nstate, TPZCompMesh *cmes { TPZMaterial *mat = cmesh->FindMaterial(matid); if(mat) return; - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(matid); + auto *nullmat = new TPZNullMaterial(matid); nullmat->SetDimension(dim); nullmat->SetNStateVariables(nstate); cmesh->InsertMaterialObject(nullmat); @@ -720,16 +720,16 @@ void TPZCreateMultiphysicsSpace::InsertLagranceMaterialObjects(TPZMultiphysicsCo { if(fSpaceType == EH1Hybrid) { - TPZLagrangeMultiplierCS<> *lag1 = new TPZLagrangeMultiplierCS<>(fH1Hybrid.fLagrangeMatid.first, fDimension-1, 1); + auto *lag1 = new TPZLagrangeMultiplier(fH1Hybrid.fLagrangeMatid.first, fDimension-1, 1); mphys->InsertMaterialObject(lag1); - TPZLagrangeMultiplierCS<> *lag2 = new TPZLagrangeMultiplierCS<>(fH1Hybrid.fLagrangeMatid.second, fDimension-1, 1); + auto *lag2 = new TPZLagrangeMultiplier(fH1Hybrid.fLagrangeMatid.second, fDimension-1, 1); lag2->SetMultiplier(-1.); mphys->InsertMaterialObject(lag2); } else if (fSpaceType == EH1HybridSquared) { - TPZLagrangeMultiplierCS<> *lag1 = new TPZLagrangeMultiplierCS<>(fH1Hybrid.fLagrangeMatid.first, fDimension-1, 1); + auto *lag1 = new TPZLagrangeMultiplier(fH1Hybrid.fLagrangeMatid.first, fDimension-1, 1); mphys->InsertMaterialObject(lag1); - TPZLagrangeMultiplierCS<> *lag2 = new TPZLagrangeMultiplierCS<>(fH1Hybrid.fSecondLagrangeMatid, fDimension-1, 1); + auto *lag2 = new TPZLagrangeMultiplier(fH1Hybrid.fSecondLagrangeMatid, fDimension-1, 1); lag2->SetMultiplier(-1.); mphys->InsertMaterialObject(lag2); } @@ -762,7 +762,7 @@ void TPZCreateMultiphysicsSpace::AddGeometricWrapElements() // first fMatWrapId TPZGeoElSide neighbour = gelside.Neighbour(); int neighmat = neighbour.Element()->MaterialId(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { if(neighmat == fH1Hybrid.fMatWrapId) { @@ -929,7 +929,7 @@ void TPZCreateMultiphysicsSpace::AssociateElements(TPZCompMesh *cmesh, TPZVec connectlist; cel->BuildConnectList(connectlist); for (auto cindex : connectlist) { -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG if (groupindex[cindex] != -1) { DebugStop(); } diff --git a/ErrorEstimation/TPZHDivErrorEstimator.cpp b/ErrorEstimation/TPZHDivErrorEstimator.cpp index 3c9d4988..2c42b8ac 100644 --- a/ErrorEstimation/TPZHDivErrorEstimator.cpp +++ b/ErrorEstimation/TPZHDivErrorEstimator.cpp @@ -7,25 +7,29 @@ // #include "TPZHDivErrorEstimator.h" +#include "DarcyFlow/TPZMixedDarcyFlow.h" +#include "TPZAnalysis.h" +#include "TPZBndCond.h" +#include "TPZCompMeshTools.h" +#include "TPZElementMatrixT.h" #include "TPZGeoElSideAncestors.h" #include "TPZGeoElSidePartition.h" #include "TPZInterfaceEl.h" -#include "TPZMixedHdivErrorEstimate.h" #include "TPZNullMaterial.h" -#include "TPZLinearAnalysis.h" +#include "TPZVTKGeoMesh.h" #include "pzbuildmultiphysicsmesh.h" #include "pzcmesh.h" #include "pzcompel.h" #include "pzcondensedcompel.h" #include "pzelementgroup.h" -#include "TPZElementMatrixT.h" +#include "pzelmat.h" #include "pzintel.h" -#include "pzsubcmesh.h" -#include "pzstepsolver.h" -#include "TPZVTKGeoMesh.h" #include "pzmultiphysicscompel.h" -#include "TPZHDivErrorEstimateMaterial.h" -#include "TPZCompMeshTools.h" +#include "pzstepsolver.h" +#include "pzsubcmesh.h" +#include +#include +#include #ifdef LOG4CXX static LoggerPtr logger(Logger::getLogger("HDivErrorEstimator")); @@ -71,22 +75,16 @@ void TPZHDivErrorEstimator::ComputeErrors(TPZVec&errorVec, TPZVec& e an.SetExact(fExact->ExactSolution()); } - int64_t nErrorCols = 6; - errorVec.resize(nErrorCols); - for (int64_t i = 0; i < nErrorCols; i++) { - errorVec[i] = 0; - } - + an.SetStep(fAdaptivityStep); + int64_t nelem = fPostProcMesh.NElements(); fPostProcMesh.LoadSolution(fPostProcMesh.Solution()); fPostProcMesh.ExpandSolution(); fPostProcMesh.ElementSolution().Redim(nelem, 5); - for(int64_t el = 0; el(cel); - if(subc) - { + auto *subc = dynamic_cast(cel); + if(subc) { int64_t nelsub = subc->NElements(); subc->ElementSolution().Redim(nelsub, 5); } @@ -94,12 +92,13 @@ void TPZHDivErrorEstimator::ComputeErrors(TPZVec&errorVec, TPZVec& e // Calculate error and store in element solution bool store_error = true; + an.SetThreadsForError(12); an.PostProcessError(errorVec, store_error); - + TPZCompMeshTools::UnCondensedElements(&fPostProcMesh); TPZCompMeshTools::UnGroupElements(&fPostProcMesh); - if (fExact) ComputeEffectivityIndices(); + ComputeEffectivityIndices(); if (!vtkPath.empty()) { PostProcessing(an, vtkPath); @@ -109,15 +108,19 @@ void TPZHDivErrorEstimator::ComputeErrors(TPZVec&errorVec, TPZVec& e for (REAL & elementerror : elementErrors) { elementerror = 0; } + + constexpr int64_t nErrorCols = 6; + errorVec.Resize(nErrorCols); + + TPZFMatrix &elsol = fPostProcMesh.ElementSolution(); for (int64_t i = 0; i < nelem; i++) { TPZCompEl *cel = fPostProcMesh.Element(i); if (!cel) continue; TPZGeoEl* gel = fPostProcMesh.Element(i)->Reference(); if (!gel) continue; - TPZFMatrix &elsol = fPostProcMesh.ElementSolution(); elementErrors[gel->Index()] = elsol(i, 3); } - + } void TPZHDivErrorEstimator::PostProcessing(TPZAnalysis &an, std::string &out) { @@ -142,7 +145,8 @@ void TPZHDivErrorEstimator::PostProcessing(TPZAnalysis &an, std::string &out) { vecnames.Push("FluxFem"); vecnames.Push("FluxReconstructed"); scalnames.Push("POrder"); - + scalnames.Push("Permeability"); + int dim = fPostProcMesh.Reference()->Dimension(); an.DefineGraphMesh(dim, scalnames, vecnames, out); @@ -193,19 +197,21 @@ TPZCompMesh *TPZHDivErrorEstimator::CreatePressureMesh() { // Insert BC materials in pressure reconstruction mesh std::set bcMatIDs = GetBCMatIDs(&fPostProcMesh); + // Insert dummy volumetric material since we need ContributeBC + auto *dummy_mat = new TPZDarcyFlow(-999, fPostProcMesh.Dimension()); + for (auto bcID : bcMatIDs) { TPZMaterial *mat = mult->FindMaterial(bcID); - TPZBndCondT *bc = dynamic_cast *>(mat); + auto *bc = dynamic_cast *>(mat); if (!bc) DebugStop(); int volumetricMatId = bc->Material()->Id(); - TPZMaterial *pressuremat = pressureMesh->FindMaterial(volumetricMatId); + TPZMaterial * volumetricMat = pressureMesh->FindMaterial(volumetricMatId); + auto *pressuremat = dynamic_cast*>(volumetricMat); if (!pressuremat) DebugStop(); - TPZMaterialT *press = dynamic_cast *>(pressuremat); - TPZBndCondT *newbc = press->CreateBC(pressuremat, bc->Id(), bc->Type(), bc->Val1(), bc->Val2()); + TPZBndCondT *newbc = dummy_mat->CreateBC(dummy_mat, bc->Id(), bc->Type(), bc->Val1(), bc->Val2()); if (bc->HasForcingFunctionBC()) { - newbc->SetForcingFunctionBC(bc->ForcingFunctionBC()); } pressureMesh->InsertMaterialObject(newbc); @@ -246,7 +252,7 @@ TPZCompMesh *TPZHDivErrorEstimator::CreatePressureMesh() { TPZGeoEl *bcGeoEl = gmesh->Element(bcElID); // We load the volumetric element reference so it shares its connects with the BC element to be created neighCel->LoadElementReference(); - TPZInterpolatedElement *intel = dynamic_cast(neighCel); + auto *intel = dynamic_cast(neighCel); if (!intel) DebugStop(); // Get order of the neighbour side of the volumetric element int iside = neighSide.Side(); @@ -260,7 +266,7 @@ TPZCompMesh *TPZHDivErrorEstimator::CreatePressureMesh() { neighCel->Reference()->ResetReference(); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream outTXT("PostProcPressureMesh.txt"); std::ofstream outVTK("PostProcPressureMesh.vtk"); @@ -279,7 +285,7 @@ TPZCompMesh *TPZHDivErrorEstimator::CreatePressureMesh() { /// hybridized) void TPZHDivErrorEstimator::CreatePostProcessingMesh() { -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream out("OriginalFlux.txt"); fOriginal->MeshVector()[0]->Print(out); @@ -292,20 +298,13 @@ void TPZHDivErrorEstimator::CreatePostProcessingMesh() { // initialize the post processing mesh fPostProcMesh.SetReference(fOriginal->Reference()); + fPostProcMesh.ApproxSpace().SetAllCreateFunctionsMultiphysicElem(); int dim = fOriginal->Dimension(); - fOriginal->CopyMaterials(fPostProcMesh); - - { - fPostProcMesh.DeleteMaterial(fHybridizer.fHDivWrapMatid); - fPostProcMesh.DeleteMaterial(fHybridizer.fInterfaceMatid.first); - fPostProcMesh.DeleteMaterial(fHybridizer.fInterfaceMatid.second); - fPostProcMesh.DeleteMaterial(fHybridizer.fLagrangeInterface); - } - // switch the material from mixed to TPZMixedHdivErrorEstimate... - SwitchMaterialObjects(); - - TPZManVector meshvec(4, 0); + InsertPostProcMaterials(); + + auto & meshvec = fPostProcMesh.MeshVector(); + meshvec.resize(4); meshvec[0] = 0; meshvec[2] = fOriginal->MeshVector()[0];//flux meshvec[3] = fOriginal->MeshVector()[1];//potential @@ -335,7 +334,7 @@ void TPZHDivErrorEstimator::CreatePostProcessingMesh() { if (dim == 3) { CreateEdgeSkeletonMesh(meshvec[1]); } -#ifdef ERRORESTIMATION_DEBUG2 +#ifdef PZDEBUG2 { if(fPostProcesswithHDiv) { @@ -347,8 +346,8 @@ void TPZHDivErrorEstimator::CreatePostProcessingMesh() { } #endif - - TPZManVector active(4, 0); + auto & active = fPostProcMesh.GetActiveApproximationSpaces(); + active.Resize(4, 0); active[1] = 1; if(fPostProcesswithHDiv) @@ -369,7 +368,7 @@ void TPZHDivErrorEstimator::CreatePostProcessingMesh() { // construction of the multiphysics mesh //cria elementos de interface fHybridizer.CreateInterfaceElements(&fPostProcMesh); - fHybridizer.GroupandCondenseElements(&fPostProcMesh); + TPZHybridizeHDiv::GroupandCondenseElements(&fPostProcMesh); fPostProcMesh.CleanUpUnconnectedNodes(); } else { @@ -378,7 +377,7 @@ void TPZHDivErrorEstimator::CreatePostProcessingMesh() { ComputePressureWeights(); -#ifdef ERRORESTIMATION_DEBUG2 +#ifdef PZDEBUG2 { std::ofstream out("multiphysicsgrouped.txt"); fPostProcMesh.Print(out); @@ -396,16 +395,16 @@ void TPZHDivErrorEstimator::ComputeElementStiffnesses() { for (auto cel:fPostProcMesh.ElementVec()) { if (!cel) continue; - TPZCondensedCompEl *condense = dynamic_cast(cel); + auto *condense = dynamic_cast(cel); if (condense) { // for Mark proposal ek correspond to local dirichlet problem condense->Assemble(); } - TPZSubCompMesh *subcmesh = dynamic_cast(cel); + auto *subcmesh = dynamic_cast(cel); if (subcmesh) { subcmesh->Assemble(); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG if(subcmesh && condense) { DebugStop(); @@ -427,7 +426,7 @@ void TPZHDivErrorEstimator::IncreaseSideOrders(TPZCompMesh *mesh) { if (gel->Dimension() != dim) { continue; } - TPZInterpolatedElement *intel = dynamic_cast (cel); + auto *intel = dynamic_cast (cel); int nc = cel->NConnects(); int order = cel->Connect(nc - 1).Order(); int nsides = gel->NSides(); @@ -483,12 +482,12 @@ void TPZHDivErrorEstimator::IncreasePressureSideOrders(TPZCompMesh *cmesh) { int maxOrder = 0; for (int ineigh = 0; ineigh < nneigh; ineigh++) { - TPZInterpolatedElement *intelS = dynamic_cast(celstack[ineigh].Element()); + auto *intelS = dynamic_cast(celstack[ineigh].Element()); int orderEl = intelS->GetPreferredOrder(); maxOrder = (orderEl > maxOrder) ? orderEl : maxOrder; } - TPZInterpolatedElement *intel = dynamic_cast (cel); + auto *intel = dynamic_cast (cel); int nsides = gel->NSides(); int ncorner = gel->NCornerNodes(); @@ -521,7 +520,7 @@ void TPZHDivErrorEstimator::CreateEdgeSkeletonMesh(TPZCompMesh *pressuremesh) { if (pressuremesh->MaterialVec().find(fPressureSkeletonMatId) != pressuremesh->MaterialVec().end()) { DebugStop(); } - TPZNullMaterial<> *nullmat = new TPZNullMaterial<>(fPressureSkeletonMatId); + auto *nullmat = new TPZNullMaterial(fPressureSkeletonMatId); pressuremesh->InsertMaterialObject(nullmat); int dim = fPostProcMesh.Dimension(); int64_t nel = pressuremesh->NElements(); @@ -530,7 +529,7 @@ void TPZHDivErrorEstimator::CreateEdgeSkeletonMesh(TPZCompMesh *pressuremesh) { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressuremesh->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = cel->Reference(); if (!gel) DebugStop(); @@ -550,7 +549,7 @@ void TPZHDivErrorEstimator::CreateEdgeSkeletonMesh(TPZCompMesh *pressuremesh) { gelpressures[createdelement->Index()] = polynomialorder; } else { int64_t gelindex = hasneigh.Element()->Index(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG if (gelpressures.find(gelindex) == gelpressures.end()) { DebugStop(); } @@ -577,7 +576,7 @@ void TPZHDivErrorEstimator::CreateEdgeSkeletonMesh(TPZCompMesh *pressuremesh) { int64_t celindex = -1; pressuremesh->SetDefaultOrder(polynomialorder); cel = pressuremesh->ApproxSpace().CreateCompEl(gel, *pressuremesh, celindex); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { TPZInterpolatedElement *intel = dynamic_cast(cel); if (!intel) DebugStop(); @@ -614,7 +613,7 @@ void TPZHDivErrorEstimator::RestrainSmallEdges(TPZCompMesh *pressuremesh) { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressuremesh->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = cel->Reference(); if (!gel) DebugStop(); @@ -627,7 +626,7 @@ void TPZHDivErrorEstimator::RestrainSmallEdges(TPZCompMesh *pressuremesh) { bool onlyinterpolated = true; TPZCompElSide large_celside = gelside.LowerLevelCompElementList2(onlyinterpolated); if (large_celside) { - TPZInterpolatedElement *largeintel = dynamic_cast(large_celside.Element()); + auto *largeintel = dynamic_cast(large_celside.Element()); if (!largeintel) DebugStop(); int largeside = large_celside.Side(); intel->RestrainSide(nsides - 1, largeintel, largeside); @@ -636,9 +635,9 @@ void TPZHDivErrorEstimator::RestrainSmallEdges(TPZCompMesh *pressuremesh) { TPZGeoElSide gelside_small(gel, side); TPZCompElSide celside_restraint = gelside_small.LowerLevelCompElementList2(onlyinterpolated); if (celside_restraint) { - TPZInterpolatedElement *largeintel = dynamic_cast(celside_restraint.Element()); + largeintel = dynamic_cast(celside_restraint.Element()); if (!largeintel) DebugStop(); - int largeside = large_celside.Side(); + largeside = large_celside.Side(); intel->RestrainSide(side, largeintel, largeside); } } @@ -659,7 +658,7 @@ void TPZHDivErrorEstimator::AdjustNeighbourPolynomialOrders(TPZCompMesh *pressur for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressureHybrid->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() >= dim) continue; @@ -671,7 +670,7 @@ void TPZHDivErrorEstimator::AdjustNeighbourPolynomialOrders(TPZCompMesh *pressur for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressureHybrid->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() >= dim) continue; @@ -703,7 +702,7 @@ void TPZHDivErrorEstimator::AdjustNeighbourPolynomialOrders(TPZCompMesh *pressur for (int ieq = 0; ieq < nequal; ieq++) { TPZCompEl *celneigh = celstack[ieq].Element(); int celside = celstack[ieq].Side(); - TPZInterpolatedElement *intelneigh = dynamic_cast(celneigh); + auto *intelneigh = dynamic_cast(celneigh); if (!intelneigh) DebugStop(); int nneighconnects = intelneigh->NSideConnects(celside); int neighporder = intelneigh->SideConnect(nneighconnects - 1, celside).Order(); @@ -718,7 +717,7 @@ void TPZHDivErrorEstimator::AdjustNeighbourPolynomialOrders(TPZCompMesh *pressur for (int ieq = 0; ieq < nequal; ieq++) { TPZCompEl *celneigh = celstack[ieq].Element(); int celside = celstack[ieq].Side(); - TPZInterpolatedElement *intelneigh = dynamic_cast(celneigh); + auto *intelneigh = dynamic_cast(celneigh); if (!intelneigh) DebugStop(); int nneighconnects = intelneigh->NSideConnects(celside); int neighporder = intelneigh->SideConnect(nneighconnects - 1, celside).Order(); @@ -743,7 +742,7 @@ void TPZHDivErrorEstimator::AdjustNeighbourPolynomialOrders(TPZCompMesh *pressur int porder = it.second; TPZCompEl *cel = pressureHybrid->Element(index); if (!cel) DebugStop(); - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); intel->SetSideOrder(side, porder); } @@ -767,7 +766,7 @@ void TPZHDivErrorEstimator::ComputeAveragePressures(int target_dim) { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressure_mesh->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() != target_dim + 1 && gel->MaterialId() != fPressureSkeletonMatId) continue; @@ -786,14 +785,14 @@ void TPZHDivErrorEstimator::ComputeAveragePressures(int target_dim) { int matid = gel->MaterialId(); TPZMaterial *mat = pressure_mesh->FindMaterial(matid); // TODO change this. Look for matIDs in bcMatIds instead. Only cast in debug mode for further checks - TPZBndCond *bc = dynamic_cast(mat); + auto *bc = dynamic_cast(mat); if (bc) continue; // Skip calculation if the element is a small skeleton bool largeSideExists = false; if (cel->Connect(0).HasDependency()) largeSideExists = true; -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG int nsides = gel->NSides(); TPZGeoElSide side(gel, nsides - 1); TPZGeoElSideAncestors ancestors(side); @@ -844,7 +843,7 @@ void TPZHDivErrorEstimator::ComputeBoundaryL2Projection(int target_dim){ int matid = gel->MaterialId(); TPZMaterial *mat = pressuremesh->FindMaterial(matid); - TPZBndCond *bc = dynamic_cast(mat); + auto *bc = dynamic_cast(mat); if (!bc || (bc->Type() != 0)) continue; cel->CalcStiff(ekbc, efbc); @@ -855,7 +854,7 @@ void TPZHDivErrorEstimator::ComputeBoundaryL2Projection(int target_dim){ TPZConnect &c = cel->Connect(ic); int64_t seqnum = c.SequenceNumber(); int64_t pos = pressuremesh->Block().Position(seqnum); - int ndof = c.NShape() * c.NState(); + uint ndof = c.NShape() * c.NState(); for (int idf = 0; idf < ndof; idf++) { mesh_sol(pos + idf, 0) = efbc.fMat(count++); } @@ -877,12 +876,12 @@ void TPZHDivErrorEstimator::ComputeAverage(TPZCompMesh *pressuremesh, int64_t ie if (!cel) DebugStop(); TPZGeoEl *gel = cel->Reference(); if (!gel) DebugStop(); - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); int nc = cel->NConnects(); int order = cel->Connect(nc - 1).Order(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG for (int ic = 0; ic < nc; ic++) { if (cel->Connect(ic).HasDependency()) DebugStop(); } @@ -921,7 +920,7 @@ void TPZHDivErrorEstimator::ComputeAverage(TPZCompMesh *pressuremesh, int64_t ie } } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG for (int ineigh = 0; ineigh < volumeNeighSides.size(); ineigh++) { if (volumeNeighSides[ineigh].Element()->Dimension() != dim) DebugStop(); } @@ -1048,7 +1047,7 @@ void TPZHDivErrorEstimator::ComputeAverage(TPZCompMesh *pressuremesh, int64_t ie TPZConnect &c = cel->Connect(ic); int64_t seqnum = c.SequenceNumber(); int64_t pos = pressuremesh->Block().Position(seqnum); - int ndof = c.NShape() * c.NState(); + uint ndof = c.NShape() * c.NState(); for (int idf = 0; idf < ndof; idf++) { mesh_sol(pos + idf, 0) = L2Rhs(count++); } @@ -1080,7 +1079,7 @@ void TPZHDivErrorEstimator::TransferEdgeSolution() { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressureHybrid->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() > 2) continue; @@ -1091,7 +1090,7 @@ void TPZHDivErrorEstimator::TransferEdgeSolution() { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressureHybrid->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = cel->Reference(); if (!gel) DebugStop(); @@ -1122,7 +1121,7 @@ void TPZHDivErrorEstimator::TransferEdgeSolution() { int nshape_edge = pressureHybrid->Block().Size(edge_seqnum); for (int ieq = 0; ieq < nequal; ieq++) { TPZCompEl *celneigh = equal[ieq].Element(); - TPZInterpolatedElement *intelneigh = dynamic_cast(celneigh); + auto *intelneigh = dynamic_cast(celneigh); TPZGeoEl *neighgel = intelneigh->Reference(); if (!intelneigh) DebugStop(); if (neighgel->MaterialId() != fHybridizer.fLagrangeInterface) { @@ -1151,7 +1150,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverages() { TPZGeoMesh *gmesh = pressure_mesh->Reference(); gmesh->ResetReference(); int dim = gmesh->Dimension(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG TPZMaterial *mat = pressure_mesh->FindMaterial(fPressureSkeletonMatId); if (!mat) DebugStop(); #endif @@ -1160,7 +1159,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverages() { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressure_mesh->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() != dim - 1) continue; @@ -1173,7 +1172,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverages() { TPZStack nodesToImposeSolution; for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressure_mesh->Element(el); - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) { continue; } @@ -1225,7 +1224,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverages() { if (neigh_celside.Reference().Dimension() != 0) continue; // Get solution of the neighbour - TPZInterpolatedElement *neigh_intel = dynamic_cast(neigh_celside.Element()); + auto *neigh_intel = dynamic_cast(neigh_celside.Element()); if (!neigh_intel) DebugStop(); int64_t neigh_conindex = neigh_intel->ConnectIndex(neigh_celside.Side()); @@ -1240,7 +1239,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverages() { } // Set solution to given connect - TPZInterpolatedElement *intel = dynamic_cast(node_celside.Element()); + auto *intel = dynamic_cast(node_celside.Element()); if (!intel) DebugStop(); int side = node_gelside.Side(); @@ -1288,7 +1287,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) // and the solution of that connect. The weight of Dirichlet condition is // higher and will be used later to impose the value of the BC in the // connects when needed - std::map>> connects; + std::map> connects; for (int elc = 0; elc < celstack.size(); elc++) { TPZCompElSide celside = celstack[elc]; TPZGeoElSide gelside = celside.Reference(); @@ -1296,14 +1295,14 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) continue; } - TPZInterpolatedElement *intel = dynamic_cast(celside.Element()); + auto *intel = dynamic_cast(celside.Element()); if (!intel) DebugStop(); int64_t index = intel->Index(); REAL weight = fPressureweights[index]; if (IsZero(weight)) { - TPZMaterial *mat = intel->Material(); - TPZBndCond *bc = dynamic_cast(mat); + mat = intel->Material(); + auto *bc = dynamic_cast(mat); if (!bc) DebugStop(); continue; @@ -1319,10 +1318,7 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) int64_t seqnum = c.SequenceNumber(); if (c.NState() != nstate || c.NShape() != 1) DebugStop(); TPZManVector pt(0), x(3); - TPZManVector sol(nstate, 0.); - for (int istate = 0; istate < nstate; istate++) { - sol[istate] = solMatrix.at(block.at(seqnum, 0, istate, 0)); - } + STATE sol = solMatrix.at(block.at(seqnum, 0, 0, 0)); #ifdef LOG4CXX if(logger->isDebugEnabled()) { @@ -1336,11 +1332,11 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) TPZManVector averageSol(nstate, 0); REAL sum_weight = 0.; - for (auto it = connects.begin(); it != connects.end(); it++) { - REAL weight = it->second.first; + for (auto & connect : connects) { + REAL weight = connect.second.first; sum_weight += weight; for (int istate = 0; istate < nstate; istate++) { - averageSol[istate] += it->second.second[istate] * weight; + averageSol[istate] += connect.second.second * weight; } } @@ -1358,8 +1354,8 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) LOGPZ_DEBUG(logger, sout.str()) } #endif - for (auto it = connects.begin(); it != connects.end(); it++) { - int64_t conindex = it->first; + for (auto & connect : connects) { + int64_t conindex = connect.first; TPZConnect &c = pressure_mesh->ConnectVec()[conindex]; int64_t seqnum = c.SequenceNumber(); @@ -1385,28 +1381,25 @@ void TPZHDivErrorEstimator::ComputeNodalAverage(TPZCompElSide &node_celside) void TPZHDivErrorEstimator::ComputeEffectivityIndices(TPZSubCompMesh *subcmesh) { TPZFMatrix &elsol = subcmesh->ElementSolution(); int64_t nrows = elsol.Rows(); - int64_t ncols = 5; // subcmesh->ElementSolution().Cols(); - - if (subcmesh->ElementSolution().Cols() != 7) { - // TODO I made some changes to be able to run the code again. - // Sometimes subcmesh->ElementSolution().Cols() equals 5 sometimes it's already 7. - // I'm not sure if this behaviour is expected. - subcmesh->ElementSolution().Resize(nrows, 7); + if (subcmesh->ElementSolution().Cols() != 5) { + // DebugStop(); } + subcmesh->ElementSolution().Resize(nrows, 7); int64_t nel = subcmesh->NElements(); TPZManVector errors(5, 0.); for (int64_t el = 0; el < nel; el++) { - for (int i = 0; i < ncols; i++) { + for (int i = 0; i < 5; i++) { errors[i] += elsol(el, i) * elsol(el, i); } } - for (int i = 0; i < ncols; i++) { + for (int i = 0; i < 5; i++) { errors[i] = sqrt(errors[i]); } + // Substitute original errors by a constant in each subcmesh for (int64_t el = 0; el < nel; el++) { - for (int i = 0; i < ncols; i++) { + for (int i = 0; i < 5; i++) { elsol(el, i) = errors[i]; } } @@ -1429,11 +1422,11 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices(TPZSubCompMesh *subcmesh) } if (abs(ErrorEstimate) < tol) { - elsol(el, ncols + i / 2) = 1.; + elsol(el, 5 + i / 2) = 1.; } else { REAL EfIndex = (ErrorEstimate + oscillatoryterm) / ErrorExact; - elsol(el, ncols + i / 2) = EfIndex; + elsol(el, 5 + i / 2) = EfIndex; } } } @@ -1441,12 +1434,13 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices(TPZSubCompMesh *subcmesh) /// compute the effectivity indices of the pressure error and flux error and store in the element solution void TPZHDivErrorEstimator::ComputeEffectivityIndices() { - /**The ElementSolution() is a matrix with 4 cols, + /**The ElementSolution() is a matrix with 5 cols, col 0: pressure exact error col 1: pressure estimate error col 2: flux exact error col 3: flux estimate error - Is increased 2 cols on ElementSolution() to store the efectivity index for pressure and flux + col 4: residual/oscillatory term + Is increased 2 cols on ElementSolution() to store the effectivity index for pressure and flux **/ TPZCompMesh *cmesh = &fPostProcMesh; @@ -1479,8 +1473,9 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices() { TPZCompEl *cel = cmesh->Element(el); if (!cel) continue; - TPZSubCompMesh *subcmesh = dynamic_cast(cel); + auto *subcmesh = dynamic_cast(cel); if (subcmesh) { + //continue; ComputeEffectivityIndices(subcmesh); } TPZGeoEl *gel = cel->Reference(); @@ -1559,9 +1554,9 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices() { TPZCompEl *cel = cmesh->Element(el); if (!cel) continue; - TPZSubCompMesh *subcmesh = dynamic_cast(cel); + auto *subcmesh = dynamic_cast(cel); if (subcmesh) { - ComputeEffectivityIndices(subcmesh); + //ComputeEffectivityIndices(subcmesh); } TPZGeoEl *gel = cel->Reference(); if (!gel) continue; @@ -1584,9 +1579,9 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices() { } #endif - TPZGeoEl *gel = cel->Reference(); + gel = cel->Reference(); - REAL hk = gel->CharacteristicSize(); + hk = gel->CharacteristicSize(); REAL oscilatorytherm = 0; if (i == 2) { @@ -1621,11 +1616,11 @@ void TPZHDivErrorEstimator::ComputeEffectivityIndices() { /// returns true if the material associated with the element is a boundary condition /// and if the boundary condition is dirichlet type -bool TPZHDivErrorEstimator::IsDirichletCondition(TPZGeoElSide gelside) { +bool TPZHDivErrorEstimator::IsDirichletCondition(const TPZGeoElSide& gelside) { TPZGeoEl *gel = gelside.Element(); int matid = gel->MaterialId(); TPZMaterial *mat = fPostProcMesh.FindMaterial(matid); - TPZBndCond *bc = dynamic_cast(mat); + auto *bc = dynamic_cast(mat); if (!bc) return false; int typ = bc->Type(); // std::cout<<"type "<< typ<<"\n"; @@ -1635,11 +1630,11 @@ bool TPZHDivErrorEstimator::IsDirichletCondition(TPZGeoElSide gelside) { void TPZHDivErrorEstimator::PotentialReconstruction() { - if (fPostProcMesh.MeshVector().size()) { + if (fPostProcMesh.MeshVector().size() && !fHasReferenceSolution) { DebugStop(); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG // Create directories to store debugging files std::string command; command = "mkdir -p ReconstructionSteps"; @@ -1650,7 +1645,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { CreatePostProcessingMesh(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { REAL presnorm = Norm(fPostProcMesh.MeshVector()[1]->Solution()); if(std::isnan(presnorm)) DebugStop(); @@ -1665,7 +1660,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { //BoundaryPressureProjection(pressuremesh, target_dim); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { PlotPressureSkeleton("ReconstructionSteps/SkelBoundaryProjection"); } @@ -1679,7 +1674,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { ComputeAveragePressures(1); } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { PlotPressureSkeleton("ReconstructionSteps/SkelInterfaceAverage"); } @@ -1687,7 +1682,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { ComputeNodalAverages(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { PlotPressureSkeleton("ReconstructionSteps/SkelNodalAverage"); } @@ -1711,7 +1706,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { std::ofstream out("ReconstructionSteps/MFMeshBeforeManualTransfer.txt"); fPostProcMesh.Print(out); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { PlotState("ReconstructionSteps/VolumeMFPressureAfterTransferFromMeshes", 2, &fPostProcMesh, false); PlotState("ReconstructionSteps/VolumePressureAfterTransferFromMeshes", 2, fPostProcMesh.MeshVector()[1]); @@ -1748,7 +1743,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG VerifySolutionConsistency(PressureMesh()); #endif @@ -1759,7 +1754,7 @@ void TPZHDivErrorEstimator::PotentialReconstruction() { } } -void TPZHDivErrorEstimator::PlotPressureSkeleton(const std::string &filename, bool reconstructed) { +[[maybe_unused]] void TPZHDivErrorEstimator::PlotPressureSkeleton(const std::string &filename, bool reconstructed) { TPZCompMesh *pressure = nullptr; @@ -1767,11 +1762,10 @@ void TPZHDivErrorEstimator::PlotPressureSkeleton(const std::string &filename, bo if (!reconstructed) { pressure = fOriginal->MeshVector()[1]; - scalnames.Push("State"); } else { pressure = PressureMesh(); - scalnames.Push("State"); } + scalnames.Push("State"); TPZLinearAnalysis an(pressure, false); @@ -1788,7 +1782,7 @@ void TPZHDivErrorEstimator::PlotPressureSkeleton(const std::string &filename, bo } } -void TPZHDivErrorEstimator::PlotInterfaceFluxes(const std::string &filename, bool reconstructed) { +[[maybe_unused]] void TPZHDivErrorEstimator::PlotInterfaceFluxes(const std::string &filename, bool reconstructed) { TPZCompMesh *flux_mesh = nullptr; if (reconstructed) { if (!fPostProcesswithHDiv) DebugStop(); @@ -1816,35 +1810,38 @@ void TPZHDivErrorEstimator::PlotInterfaceFluxes(const std::string &filename, boo } /// switch material object from mixed poisson to TPZMixedHdivErrorEstimate -void TPZHDivErrorEstimator::SwitchMaterialObjects() { +void TPZHDivErrorEstimator::InsertPostProcMaterials() { - for (auto mat : fPostProcMesh.MaterialVec()) { - TPZMixedPoisson *mixpoisson = dynamic_cast(mat.second); - if (mixpoisson) { - TPZMixedPoisson *newmat; + for (auto it : fOriginal->MaterialVec()) { + auto *orig_mat = dynamic_cast(it.second); + if (orig_mat) { + TPZMixedDarcyFlow *new_mat; if (fPostProcesswithHDiv) { - newmat = new TPZMixedHDivErrorEstimate(*mixpoisson); + new_mat = new TPZMixedHDivErrorEstimate(*orig_mat); } else { - newmat = new TPZHDivErrorEstimateMaterial(*mixpoisson); + auto ee_mat = new TPZHDivErrorEstimateMaterial(*orig_mat); + ee_mat->SetReferenceSolution(fHasReferenceSolution); + new_mat = ee_mat; } - - if (mixpoisson->HasForcingFunction()) { - newmat->SetForcingFunction(mixpoisson->ForcingFunction(), - mixpoisson->ForcingFunctionPOrder()); + if (orig_mat->HasForcingFunction()) { + new_mat->SetForcingFunction(orig_mat->ForcingFunction(), orig_mat->ForcingFunctionPOrder()); } - if (mixpoisson->HasExactSol()) { - newmat->SetExactSol(mixpoisson->ExactSol(), - mixpoisson->PolynomialOrderExact()); + if (orig_mat->HasExactSol()) { + new_mat->SetExactSol(orig_mat->ExactSol(), orig_mat->PolynomialOrderExact()); } + fPostProcMesh.InsertMaterialObject(new_mat); + } + } - for (auto bcmat : fPostProcMesh.MaterialVec()) { - TPZBndCond *bc = dynamic_cast(bcmat.second); - if (bc) { - bc->SetMaterial(newmat); - } - } - fPostProcMesh.MaterialVec()[newmat->Id()] = newmat; - delete mixpoisson; + for (auto it : fOriginal->MaterialVec()) { + auto *orig_bc = dynamic_cast(it.second); + if (orig_bc) { + it.second->Clone(fPostProcMesh.MaterialVec()); + auto *new_mat = fPostProcMesh.FindMaterial(orig_bc->Material()->Id()); + if (!new_mat) DebugStop(); + auto *new_bc = dynamic_cast(fPostProcMesh.FindMaterial(orig_bc->Id())); + if (!new_bc) DebugStop(); + new_bc->SetMaterial(new_mat); } } } @@ -1954,7 +1951,7 @@ void TPZHDivErrorEstimator::VerifySolutionConsistency(TPZCompMesh *cmesh) { #endif // Checks pressure value on these nodes - TPZInterpolatedElement *intel = dynamic_cast(cneighbour.Element()); + auto *intel = dynamic_cast(cneighbour.Element()); if (!intel) DebugStop(); } } @@ -1987,17 +1984,17 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { TPZGeoElSide skelSide(gel, gel->NSides() - 1); TPZStack compNeighSides; skelSide.EqualLevelCompElementList(compNeighSides, 1, 0); - if (compNeighSides.size() == 1) { - TPZCompElSide large = skelSide.LowerLevelCompElementList2(true); - if (large) { - std::cout << "Gel: " << gel->Index() << " Side: " << skelSide.Side() << '\n'; - //compNeighSides.Push(large); - } - } - + //if (compNeighSides.size() == 1) { + // TPZCompElSide large = skelSide.LowerLevelCompElementList2(true); + // if (large) { + // std::cout << "Gel: " << gel->Index() << " Side: " << skelSide.Side() << '\n'; + // //compNeighSides.Push(large); + // } + //} + for (int i = 0; i < compNeighSides.size(); i++) { TPZCompEl *neighCel = compNeighSides[i].Element(); - TPZInterpolatedElement *neighIntEl = dynamic_cast(neighCel); + auto *neighIntEl = dynamic_cast(neighCel); if (!neighIntEl) DebugStop(); int sideNum = compNeighSides[i].Side(); @@ -2021,7 +2018,7 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { if (!gel) DebugStop(); if (gel->Dimension() != fPostProcMesh.Dimension()) continue; - TPZMultiphysicsElement *mcel = dynamic_cast(cel); + auto *mcel = dynamic_cast(cel); if (!mcel) DebugStop(); std::set elementsToGroup; @@ -2036,7 +2033,7 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { TPZCompElSide compSide = compNeighSides[i]; TPZCompEl *neighCel = compSide.Element(); TPZMaterial *mat = neighCel->Material(); - TPZBndCond *bc = dynamic_cast(mat); + auto *bc = dynamic_cast(mat); if (bc) { elementsToGroup.insert(mcel->Index()); elementsToGroup.insert(neighCel->Index()); @@ -2044,9 +2041,9 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { } } - if (elementsToGroup.size()) { + if (!elementsToGroup.empty()) { int64_t index; - TPZElementGroup *elGroup = new TPZElementGroup(fPostProcMesh, index); + auto *elGroup = new TPZElementGroup(fPostProcMesh, index); for (const auto &it : elementsToGroup) { elGroup->AddElement(fPostProcMesh.Element(it)); } @@ -2061,7 +2058,7 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { } } -#ifdef ERRORESTIMATION_DEBUG2 +#ifdef PZDEBUG2 { std::ofstream txtPostProcMesh("PostProcMeshAfterIncrementingConnects.txt"); fPostProcMesh.Print(txtPostProcMesh); @@ -2075,17 +2072,17 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { if (!cel) continue; TPZGeoEl *gel = cel->Reference(); if (!gel) { - TPZElementGroup *group = dynamic_cast(cel); + auto *group = dynamic_cast(cel); if (!group) DebugStop(); } if (gel && gel->Dimension() != fPostProcMesh.Dimension()) continue; - TPZCondensedCompEl *condense = new TPZCondensedCompEl(cel, false); + auto *condense = new TPZCondensedCompEl(cel, false); } // @TODO what is the meaning of this? phil for (auto matit : fPostProcMesh.MaterialVec()) { TPZMaterial *mat = matit.second; - TPZHDivErrorEstimateMaterial *errormat = dynamic_cast(mat); + auto *errormat = dynamic_cast(mat); if (errormat) { errormat->fNeumannLocalProblem = false; } @@ -2093,7 +2090,7 @@ void TPZHDivErrorEstimator::PrepareElementsForH1Reconstruction() { fPostProcMesh.CleanUpUnconnectedNodes(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG std::ofstream outTXT("PostProcessMeshAfterPreparingElements.txt"); fPostProcMesh.Print(outTXT); #endif @@ -2112,7 +2109,7 @@ void TPZHDivErrorEstimator::CopySolutionFromSkeleton() { for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = pressuremesh->Element(el); if (!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl *gel = intel->Reference(); if (gel->Dimension() != dim) continue; @@ -2122,7 +2119,7 @@ void TPZHDivErrorEstimator::CopySolutionFromSkeleton() { TPZGeoElSide gelside(gel, is); TPZConnect &c = intel->Connect(is); int64_t c_seqnum = c.SequenceNumber(); - int c_blocksize = c.NShape() * c.NState(); + uint c_blocksize = c.NShape() * c.NState(); TPZStack celstack; gelside.EqualLevelCompElementList(celstack, 1, 0); @@ -2131,11 +2128,11 @@ void TPZHDivErrorEstimator::CopySolutionFromSkeleton() { TPZCompElSide cneigh = celstack[ist]; TPZGeoElSide gneigh = cneigh.Reference(); if (gneigh.Element()->MaterialId() == this->fPressureSkeletonMatId || IsDirichletCondition(gneigh)) { - TPZInterpolatedElement *intelneigh = dynamic_cast(cneigh.Element()); + auto *intelneigh = dynamic_cast(cneigh.Element()); if (!intelneigh) DebugStop(); TPZConnect &con_neigh = intelneigh->Connect(cneigh.Side()); int64_t con_seqnum = con_neigh.SequenceNumber(); - int con_size = con_neigh.NState() * con_neigh.NShape(); + uint con_size = con_neigh.NState() * con_neigh.NShape(); if (con_size != c_blocksize) DebugStop(); for (int ibl = 0; ibl < con_size; ibl++) { @@ -2180,7 +2177,7 @@ void TPZHDivErrorEstimator::ComputePressureWeights() { } if (!mat) DebugStop(); - TPZBndCondT *bcmat = dynamic_cast *>(mat); + auto *bcmat = dynamic_cast *>(mat); if (bcmat) { switch(bcmat->Type()) { case 0: // Dirichlet BC @@ -2199,16 +2196,15 @@ void TPZHDivErrorEstimator::ComputePressureWeights() { DebugStop(); } } else { - TPZMixedPoisson *mixpoisson = dynamic_cast(mat); + auto *mixpoisson = dynamic_cast(mat); if (!mixpoisson) DebugStop(); - REAL perm; TPZVec xi(gel->Dimension(), 0.); gel->CenterPoint(gel->NSides() - 1, xi); TPZVec x(3, 0.); gel->X(xi, x); - perm = mixpoisson->GetPermeability(x); + STATE perm = mixpoisson->GetPermeability(x); if (IsZero(perm)) DebugStop(); this->fPressureweights[el] = perm; fMatid_weights[matid] = perm; @@ -2217,7 +2213,7 @@ void TPZHDivErrorEstimator::ComputePressureWeights() { std::cout << "Finished computing pressure weights\n"; } -void TPZHDivErrorEstimator::PlotState(const std::string& filename, int targetDim, TPZCompMesh* cmesh, bool atomic) { +[[maybe_unused]] void TPZHDivErrorEstimator::PlotState(const std::string& filename, int targetDim, TPZCompMesh* cmesh, bool atomic) { std::ofstream outTXT("PressuretoStateGraph.txt"); cmesh->Print(outTXT); @@ -2238,7 +2234,7 @@ void TPZHDivErrorEstimator::PlotState(const std::string& filename, int targetDim plotname = out.str(); } an.DefineGraphMesh(targetDim, scalnames, vecnames, plotname); - an.PostProcess(4, targetDim); + an.PostProcess(0, targetDim); } } @@ -2250,7 +2246,7 @@ void TPZHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh *pressure_mesh) { gmesh->ResetReference(); cmesh->LoadReferences(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream fileVTK("GeoMeshBeforePressureSkeleton.vtk"); TPZVTKGeoMesh::PrintGMeshVTK(gmesh, fileVTK); @@ -2294,7 +2290,7 @@ void TPZHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh *pressure_mesh) { } } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream fileVTK("GeoMeshAfterPressureSkeleton.vtk"); TPZVTKGeoMesh::PrintGMeshVTK(gmesh, fileVTK); @@ -2304,7 +2300,7 @@ void TPZHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh *pressure_mesh) { #endif // Create skeleton elements in pressure mesh - TPZNullMaterial<> *skeletonMat = new TPZNullMaterial<>(fPressureSkeletonMatId); + auto *skeletonMat = new TPZNullMaterial(fPressureSkeletonMatId); skeletonMat->SetDimension(dim - 1); pressure_mesh->InsertMaterialObject(skeletonMat); @@ -2323,13 +2319,13 @@ void TPZHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh *pressure_mesh) { } -void TPZHDivErrorEstimator::RestrainSkeletonSides(TPZCompMesh *pressure_mesh) { +void TPZHDivErrorEstimator::RestrainSkeletonSides(TPZCompMesh *pressure_mesh) const { TPZGeoMesh *gmesh = pressure_mesh->Reference(); gmesh->ResetReference(); pressure_mesh->LoadReferences(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream out("MeshBeforeRestrainSkeleton.txt"); pressure_mesh->Print(out); @@ -2352,21 +2348,18 @@ void TPZHDivErrorEstimator::RestrainSkeletonSides(TPZCompMesh *pressure_mesh) { TPZCompEl *smallCel = small.Element()->Reference(); if (!smallCel) DebugStop(); - TPZInterpolatedElement *smallIntel = dynamic_cast(smallCel); + auto *smallIntel = dynamic_cast(smallCel); if (!smallIntel) DebugStop(); TPZCompEl *largeCel = largerNeigh.Element()->Reference(); if (!largeCel) DebugStop(); - TPZInterpolatedElement *largeIntel = dynamic_cast(largeCel); + auto *largeIntel = dynamic_cast(largeCel); if (!largeIntel) DebugStop(); TPZManVector xicenter(small.Element()->Dimension(), 0.); gel->CenterPoint(nsides - 1, xicenter); TPZManVector xcenter(3, 0.); gel->X(xicenter, xcenter); - std::cout << "Restriction @ [" << xcenter << "]:" - << " Small El: " << small.Element()->Index() << ", Side: " << small.Side() - << " Large El: " << largerNeigh.Element()->Index() << ", Side: " << largerNeigh.Side() << "\n"; smallIntel->RestrainSide(small.Side(), largeIntel, largerNeigh.Side()); // Restrain subsides @@ -2378,16 +2371,13 @@ void TPZHDivErrorEstimator::RestrainSkeletonSides(TPZCompMesh *pressure_mesh) { if (subLargerNeigh.Element() == subsmall.Element()) DebugStop(); gel->CenterPoint(iside, xicenter); gel->X(xicenter, xcenter); - std::cout << "SubRestriction @ [" << xcenter << "]:" - << " Small El: " << small.Element()->Index() << ", Side: " << subsmall.Side() - << " Large El: " << largerNeigh.Element()->Index() << ", Side: " << subLargerNeigh.Side() << "\n"; smallIntel->RestrainSide(subsmall.Side(), largeIntel, subLargerNeigh.Side()); } } pressure_mesh->CleanUpUnconnectedNodes(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream out("MeshAfterRestrainSkeleton.txt"); pressure_mesh->Print(out); @@ -2424,7 +2414,7 @@ bool TPZHDivErrorEstimator::IsAdjacentToHangingNode(const TPZCompElSide &celside if (neigh_gelside.Element()->Dimension() != dim - 1) { continue; } - TPZInterpolatedElement *intel = dynamic_cast(neigh_celside.Element()); + auto *intel = dynamic_cast(neigh_celside.Element()); if (!intel) DebugStop(); int64_t conindex = intel->ConnectIndex(neigh_celside.Side()); @@ -2461,3 +2451,31 @@ std::set TPZHDivErrorEstimator::GetBCMatIDs(const TPZCompMesh* cmesh) { } return bc_mat_ids; } + +void TPZHDivErrorEstimator::SetReferenceSolution(bool hasRefSol) { + fHasReferenceSolution = hasRefSol; + + auto mat_vec = fPostProcMesh.MaterialVec(); + for (auto &mat : mat_vec) { + auto *material = mat.second; + auto *error_mat = dynamic_cast(material); + if (error_mat) { + error_mat->SetReferenceSolution(true); + } + } +} + +void TPZHDivErrorEstimator::SetReferenceSolutionMeshes(TPZCompMesh* ref_pressure, TPZCompMesh* ref_flux) { + fHasReferenceSolution = true; + auto &mesh_vec = fPostProcMesh.MeshVector(); + auto &active_spaces = fPostProcMesh.GetActiveApproximationSpaces(); + + mesh_vec.Resize(6, nullptr); + active_spaces.Resize(6, 0); + + mesh_vec[4] = ref_pressure; + mesh_vec[5] = ref_flux; + + active_spaces[4] = 0; + active_spaces[5] = 0; +} diff --git a/ErrorEstimation/TPZHDivErrorEstimator.h b/ErrorEstimation/TPZHDivErrorEstimator.h index 870f28fb..66f29750 100644 --- a/ErrorEstimation/TPZHDivErrorEstimator.h +++ b/ErrorEstimation/TPZHDivErrorEstimator.h @@ -33,8 +33,6 @@ class TPZHDivErrorEstimator { TPZVec fPressureweights; /// weights associated with material ids std::map fMatid_weights; - // object to assist in creating a hybridized version of the computational mesh - TPZHybridizeHDiv fHybridizer; // material id of the dim-1 skeleton elements int fPressureSkeletonMatId = 0; @@ -43,6 +41,14 @@ class TPZHDivErrorEstimator { /// whether the post processing mesh will be H(div) or H1 bool fPostProcesswithHDiv = false; + int fAdaptivityStep{0}; + + bool fHasReferenceSolution{false}; + +private: + // object to assist in creating a hybridized version of the computational mesh + TPZHybridizeHDiv fHybridizer; + public: explicit TPZHDivErrorEstimator(TPZMultiphysicsCompMesh &originalMesh, bool postProcWithHDiv = false) : fOriginal(&originalMesh), fPostProcMesh(nullptr), fExact(nullptr) { @@ -57,12 +63,14 @@ class TPZHDivErrorEstimator { ~TPZHDivErrorEstimator(); - //void SetProblemConfig(const ProblemConfig &cfg) { fProblemConfig = cfg; } - /// Set the analytic solution object void SetAnalyticSolution(TPZAnalyticSolution &exact) { fExact = &exact; } - void SetHybridizer(TPZHybridizeHDiv &hybridizer) { fHybridizer = hybridizer; } + [[maybe_unused]] void SetHybridizer(TPZHybridizeHDiv &hybridizer) { fHybridizer = hybridizer; } + + void SetAdaptivityStep(int step) { fAdaptivityStep = step; }; + + [[nodiscard]] int AdaptivityStep() const { return fAdaptivityStep; }; /// compute the element errors comparing the reconstructed solution based on average pressures /// with the original solution @@ -74,17 +82,21 @@ class TPZHDivErrorEstimator { /// create graphical output of estimated and true errors using the analysis void PostProcessing(TPZAnalysis &an, std::string &out); - void PlotPressureSkeleton(const std::string &filename, bool reconstructed = true); - void PlotInterfaceFluxes(const std::string &filename, bool reconstructed = true); + [[maybe_unused]] void PlotPressureSkeleton(const std::string &filename, bool reconstructed = true); + [[maybe_unused]] void PlotInterfaceFluxes(const std::string &filename, bool reconstructed = true); // Plots State solution of elements of target dimension - static void PlotState(const std::string& filename, int targetDim, TPZCompMesh* cmesh, bool atomic = true); + [[maybe_unused]] static void PlotState(const std::string& filename, int targetDim, TPZCompMesh* cmesh, bool atomic = true); - int PressureSkeletonMatId() const { return fPressureSkeletonMatId; } + [[nodiscard]] int PressureSkeletonMatId() const { return fPressureSkeletonMatId; } TPZMultiphysicsCompMesh *PostProcMesh() { return &fPostProcMesh; } TPZGeoMesh *GMesh() { return fOriginal->Reference(); } + void SetReferenceSolution(bool hasRefSol); + + void SetReferenceSolutionMeshes(TPZCompMesh* ref_pressure, TPZCompMesh* ref_flux); + protected: /// compute the pressure weights and material weights @@ -137,7 +149,7 @@ class TPZHDivErrorEstimator { /// adjust the interpolation orders so as to create an H1/2 boundary mesh // this method is called by the CreateEdgeSkeletonMesh method - void AdjustNeighbourPolynomialOrders(TPZCompMesh *pressure_mesh); + static void AdjustNeighbourPolynomialOrders(TPZCompMesh *pressure_mesh); /// restrain the edge elements that have larger elements as neighbours void RestrainSmallEdges(TPZCompMesh *pressuremesh); @@ -147,27 +159,25 @@ class TPZHDivErrorEstimator { /// compute the nodal average of all elements that share a point void ComputeNodalAverage(TPZCompElSide &node_celside); - //compute the global efectivity index using the CharacteristicSize() of element - void GlobalEffectivityIndex(); /// copy the solution from the neighbouring skeleton elements // this is a placeholder for the derived class TPZHDivErrorEstimatorH1 virtual void CopySolutionFromSkeleton(); /// switch material object from mixed poisson to TPZMixedHdivErrorEstimate - virtual void SwitchMaterialObjects(); + virtual void InsertPostProcMaterials(); /// compute the effectivity indices of the pressure error and flux error and store in the element solution void ComputeEffectivityIndices(); /// compute the effectivity indices of the pressure error and flux error and store in the element solution - void ComputeEffectivityIndices(TPZSubCompMesh *subcmesh); + static void ComputeEffectivityIndices(TPZSubCompMesh *subcmesh); /// returns true if the material associated with the element is a boundary condition /// and if the boundary condition is dirichlet type - bool IsDirichletCondition(TPZGeoElSide gelside); + bool IsDirichletCondition(const TPZGeoElSide& gelside); - void RestrainSkeletonSides(TPZCompMesh *pressure_mesh); + void RestrainSkeletonSides(TPZCompMesh *pressure_mesh) const; // Checks if the solution is in fact continuous virtual void VerifySolutionConsistency(TPZCompMesh* cmesh); @@ -180,6 +190,7 @@ class TPZHDivErrorEstimator { bool IsAdjacentToHangingNode(const TPZCompElSide &celside); static std::set GetBCMatIDs(const TPZCompMesh* cmesh); + }; #endif /* TPZHybridHDivErrorEstimator_hpp */ diff --git a/ErrorEstimation/TPZHybridH1ErrorEstimator.cpp b/ErrorEstimation/TPZHybridH1ErrorEstimator.cpp index d542a659..1395a846 100644 --- a/ErrorEstimation/TPZHybridH1ErrorEstimator.cpp +++ b/ErrorEstimation/TPZHybridH1ErrorEstimator.cpp @@ -1187,9 +1187,6 @@ void TPZHybridH1ErrorEstimator::RestrainSkeletonSides(TPZCompMesh *pressure_mesh if (subLargerNeigh.Element() == subsmall.Element()) DebugStop(); gel->CenterPoint(iside, xicenter); gel->X(xicenter, xcenter); - //std::cout << "SubRestriction @ [" << xcenter << "]:" - // << " Small El: " << small.Element()->Index() << ", Side: " << subsmall.Side() - // << " Large El: " << largerNeigh.Element()->Index() << ", Side: " << subLargerNeigh.Side() << "\n"; smallIntel->RestrainSide(subsmall.Side(), largeIntel, subLargerNeigh.Side()); } } diff --git a/ErrorEstimation/TPZMHMHDivErrorEstimator.cpp b/ErrorEstimation/TPZMHMHDivErrorEstimator.cpp index 227b390c..ca204593 100644 --- a/ErrorEstimation/TPZMHMHDivErrorEstimator.cpp +++ b/ErrorEstimation/TPZMHMHDivErrorEstimator.cpp @@ -1,15 +1,17 @@ #include "TPZMHMHDivErrorEstimator.h" +#include "TPZBndCond.h" +#include "TPZLagrangeMultiplierCS.h" #include "TPZNullMaterial.h" -#include "TPZNullMaterialCS.h" #include "TPZVTKGeoMesh.h" #include "pzintel.h" #include "pzsubcmesh.h" -#include +#include +#include #include #include -#include #include +#include #ifdef LOG4CXX static LoggerPtr logger(Logger::getLogger("HDivErrorEstimator")); @@ -21,12 +23,17 @@ void TPZMHMHDivErrorEstimator::CreatePostProcessingMesh() // initialize the post processing mesh fPostProcMesh.SetReference(fOriginal->Reference()); - fOriginal->CopyMaterials(fPostProcMesh); + //fOriginal->CopyMaterials(fPostProcMesh); // Switch the material from mixed to TPZMHMHDivErrorEstimationMaterial - SwitchMaterialObjects(); + InsertPostProcMaterials(); + + auto & meshvec = fPostProcMesh.MeshVector(); + auto & active = fPostProcMesh.GetActiveApproximationSpaces(); + if (!fHasReferenceSolution){ + meshvec.Resize(4, nullptr); + active.Resize(4, 0); + } - TPZManVector meshvec(4); - TPZManVector active(4,0); active[1] = 1; meshvec[0] = 0; @@ -46,12 +53,13 @@ void TPZMHMHDivErrorEstimator::CreatePostProcessingMesh() if (fPostProcesswithHDiv) { CreateFluxSkeletonElements(meshvec[0]); // TODO move this to InsertMaterials method - TPZNullMaterialCS<> *skeletonMat = new TPZNullMaterialCS<>(fPressureSkeletonMatId); - skeletonMat->SetDimension(GMesh()->Dimension() - 1); + auto *skeletonMat = new TPZLagrangeMultiplierCS( + fPressureSkeletonMatId, GMesh()->Dimension() - 1, 1); fPostProcMesh.InsertMaterialObject(skeletonMat); - TPZNullMaterialCS<> *wrapMat = new TPZNullMaterialCS<>(fHDivWrapMatId); - skeletonMat->SetDimension(GMesh()->Dimension() - 1); + auto *wrapMat = new TPZLagrangeMultiplierCS(fHDivWrapMatId, + GMesh()->Dimension() - + 1, 1); fPostProcMesh.InsertMaterialObject(wrapMat); } @@ -62,6 +70,7 @@ void TPZMHMHDivErrorEstimator::CreatePostProcessingMesh() } RemoveMaterialObjects(fPostProcMesh.MaterialVec()); + fPostProcMesh.ApproxSpace().SetAllCreateFunctionsMultiphysicElem(); fPostProcMesh.BuildMultiphysicsSpace(active, meshvec); if(fPostProcesswithHDiv) { @@ -73,7 +82,7 @@ void TPZMHMHDivErrorEstimator::CreatePostProcessingMesh() for (int i = 0; i < fluxMesh->NElements(); i++) { TPZCompEl *cel = fluxMesh->Element(i); if (!cel) continue; - TPZInterpolatedElement * intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) continue; TPZGeoEl * gel = cel->Reference(); @@ -123,7 +132,7 @@ void TPZMHMHDivErrorEstimator::SubStructurePostProcessingMesh() TPZCompEl *cel = gel->Reference(); if(!cel) continue; TPZCompMesh *mesh = cel->Mesh(); - TPZSubCompMesh *submesh = dynamic_cast(mesh); + auto *submesh = dynamic_cast(mesh); if (submesh) { orig_submesh_per_gel[iel] = submesh; } @@ -150,7 +159,7 @@ void TPZMHMHDivErrorEstimator::SubStructurePostProcessingMesh() auto iter = submeshmap.find(submesh_orig); if (iter == submeshmap.end()) { int64_t index; - TPZSubCompMesh *new_submesh = new TPZSubCompMesh(fPostProcMesh,index); + auto *new_submesh = new TPZSubCompMesh(fPostProcMesh,index); submeshmap[submesh_orig] = new_submesh; new_submesh_per_cel[el] = new_submesh; } else { @@ -161,7 +170,7 @@ void TPZMHMHDivErrorEstimator::SubStructurePostProcessingMesh() fPostProcMesh.ComputeNodElCon(); // TODO refactor here std::cout << "Transferring volumetric elements...\n"; - if(1) { + if(true) { TPZVec connectToSubcmesh(fPostProcMesh.NConnects(), nullptr); // transfer the elements in the submesh indicated by elementgroup @@ -212,7 +221,7 @@ void TPZMHMHDivErrorEstimator::SubStructurePostProcessingMesh() } } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG if (celDomain.size() > 1) DebugStop(); #endif @@ -300,10 +309,10 @@ void TPZMHMHDivErrorEstimator::SubStructurePostProcessingMesh() // set an analysis type for the submeshes { - int64_t nel = fPostProcMesh.NElements(); + nel = fPostProcMesh.NElements(); for (int64_t el = 0; el < nel; el++) { TPZCompEl *cel = fPostProcMesh.Element(el); - TPZSubCompMesh *sub = dynamic_cast(cel); + auto *sub = dynamic_cast(cel); if (sub) { if (fPostProcesswithHDiv) { @@ -335,7 +344,7 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateFluxMesh() for (int64_t el = 0; elElement(el); if(!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); TPZGeoEl *gel = cel->Reference(); bool isbcmat = fMHM->fMaterialBCIds.find(gel->MaterialId()) != fMHM->fMaterialBCIds.end(); // if the element is of lower dimension and is not a boundary @@ -391,7 +400,7 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateDiscontinuousPressureMesh() for (int64_t el = 0; elElement(el); if(!cel) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if(!intel) DebugStop(); TPZGeoEl *gel = cel->Reference(); int matid = gel->MaterialId(); @@ -447,7 +456,7 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateInternallyContinuousPressureMesh() if (neighbour.Element()->Dimension() != dim) { TPZMaterial *neighMat = fPostProcMesh.FindMaterial(neighbour.Element()->MaterialId()); if (neighMat) { - TPZBndCond *bc = dynamic_cast(neighMat); + auto *bc = dynamic_cast(neighMat); if (bc) { int64_t bcId = neighbour.Element()->Index(); MHMOfEachGeoEl[bcId] = {geoToMHM[el], bcId, nullptr}; @@ -462,7 +471,7 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateInternallyContinuousPressureMesh() std::sort(&MHMOfEachGeoEl[0], &MHMOfEachGeoEl[nel - 1] + 1); // Create pressure mesh - TPZCompMesh *reconstruction_pressure = new TPZCompMesh(gmesh); + auto *reconstruction_pressure = new TPZCompMesh(gmesh); // Copies volume materials original_pressure->CopyMaterials(*reconstruction_pressure); @@ -470,16 +479,24 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateInternallyContinuousPressureMesh() // Copies BC materials std::set bcMatIDs; + + auto * dummy_mat = new TPZDarcyFlow(-999, fPostProcMesh.Dimension()); + reconstruction_pressure->InsertMaterialObject(dummy_mat); for (auto it : fOriginal->MaterialVec()) { TPZMaterial *mat = it.second; - TPZBndCondT *bc = dynamic_cast *>(mat); + auto *bc = dynamic_cast *>(mat); if (bc) { int bcID = bc->Material()->Id(); - TPZMaterial *pressure_abs = original_pressure->FindMaterial(bcID); - TPZMaterialT *pressure_mat = dynamic_cast*>(pressure_abs); - TPZBndCondT *bc_mat = pressure_mat->CreateBC(pressure_mat, mat->Id(), bc->Type(), bc->Val1(), bc->Val2()); + TPZMaterial * press = original_pressure->FindMaterial(bcID); + auto *pressure_mat = dynamic_cast*>(original_pressure->FindMaterial(bcID)); + TPZBndCondT *bc_mat = dummy_mat->CreateBC(dummy_mat, mat->Id(), bc->Type(), bc->Val1(), bc->Val2()); if (fExact) { - bc_mat->SetForcingFunctionBC(fExact->ExactSolution()); + auto ff_bc = [this](const TPZVec &loc, TPZVec &rhsVal, + TPZFMatrix &matVal) { + this->fExact->Exact()->Execute(loc, rhsVal, matVal); + }; + + bc_mat->SetForcingFunctionBC(ff_bc); } reconstruction_pressure->InsertMaterialObject(bc_mat); } @@ -515,7 +532,7 @@ TPZCompMesh *TPZMHMHDivErrorEstimator::CreateInternallyContinuousPressureMesh() int64_t index; TPZCompEl *new_cel = reconstruction_pressure->CreateCompEl(gel, index); - TPZInterpolatedElement *new_intel = dynamic_cast(new_cel); + auto *new_intel = dynamic_cast(new_cel); TPZCompEl * orig_cel = std::get<2>(MHMOfEachGeoEl[i]); if (orig_cel) { @@ -597,7 +614,7 @@ void TPZMHMHDivErrorEstimator::TransferEmbeddedElements() // connect_submesh : pointer to the submesh if there is only one subdomain that owns the connect for (int64_t el=0; el(cel); + auto *sub = dynamic_cast(cel); if(!sub) continue; int nconnect = sub->NConnects(); for(int ic=0; ic(cel); + auto *sub = dynamic_cast(cel); if(sub) continue; std::set submeshes; int nconnect = cel->NConnects(); @@ -629,7 +646,7 @@ void TPZMHMHDivErrorEstimator::TransferEmbeddedElements() } if(submeshes.size() == 1) { - TPZSubCompMesh *sub = *submeshes.begin(); + sub = *submeshes.begin(); if(sub) { sub->TransferElement(&fPostProcMesh, el); @@ -639,7 +656,7 @@ void TPZMHMHDivErrorEstimator::TransferEmbeddedElements() fPostProcMesh.ComputeNodElCon(); for (int64_t el=0; el(cel); + auto *sub = dynamic_cast(cel); if(!sub) continue; sub->MakeAllInternal(); } @@ -679,7 +696,7 @@ void TPZMHMHDivErrorEstimator::ComputeAveragePressures(int target_dim) if(!gel) continue; if(gel->Dimension() == target_dim) { - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if(!intel) DebugStop(); int matId = gel->MaterialId(); @@ -718,7 +735,7 @@ void TPZMHMHDivErrorEstimator::ComputeNodalAverages() if(!cel) continue; TPZGeoEl *gel = cel->Reference(); if(!gel || !gel->Reference()) continue; - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if(!intel) DebugStop(); if (gel->Dimension() == dim-1) { int ncorner = gel->NCornerNodes(); @@ -739,7 +756,7 @@ void TPZMHMHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh * pressure_mes gmesh->ResetReference(); cmesh->LoadReferences(); -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream fileVTK("GeoMeshBeforePressureSkeleton.vtk"); TPZVTKGeoMesh::PrintGMeshVTK(gmesh, fileVTK); @@ -788,7 +805,7 @@ void TPZMHMHDivErrorEstimator::CreateSkeletonElements(TPZCompMesh * pressure_mes } } -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream fileVTK("GeoMeshAfterPressureSkeleton.vtk"); TPZVTKGeoMesh::PrintGMeshVTK(gmesh, fileVTK); @@ -802,7 +819,7 @@ void TPZMHMHDivErrorEstimator::CreateSkeletonApproximationSpace(TPZCompMesh *pre int dim = gmesh->Dimension(); // Create skeleton elements in pressure mesh - TPZNullMaterial<> *skeletonMat = new TPZNullMaterial<>(fPressureSkeletonMatId); + auto *skeletonMat = new TPZNullMaterial(fPressureSkeletonMatId); skeletonMat->SetDimension(dim - 1); pressure_mesh->InsertMaterialObject(skeletonMat); @@ -826,7 +843,7 @@ void TPZMHMHDivErrorEstimator::CopySolutionFromSkeleton() { int64_t nel = pressuremesh->NElements(); for (int64_t el = 0; el < nel; el++) { TPZCompEl* cel = pressuremesh->Element(el); - TPZInterpolatedElement *intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!cel) continue; if(!intel) DebugStop(); // load just d dimensional elements @@ -841,7 +858,7 @@ void TPZMHMHDivErrorEstimator::CopySolutionFromSkeleton() { for (int64_t el = 0; el < nel; el++) { TPZCompEl* cel = pressuremesh->Element(el); if (!cel) continue; - TPZInterpolatedElement* intel = dynamic_cast(cel); + auto *intel = dynamic_cast(cel); if (!intel) DebugStop(); TPZGeoEl* gel = cel->Reference(); // filters just (d-1) dimensional elements @@ -854,7 +871,7 @@ void TPZMHMHDivErrorEstimator::CopySolutionFromSkeleton() { TPZConnect &c = intel->Connect(is); int64_t c_gelSide_seqnum = c.SequenceNumber(); - int c_blocksize = c.NShape() * c.NState(); + uint c_blocksize = c.NShape() * c.NState(); TPZStack celstack; gelside.EqualLevelCompElementList(celstack, 1, 0); @@ -865,11 +882,11 @@ void TPZMHMHDivErrorEstimator::CopySolutionFromSkeleton() { TPZCompElSide cneigh = celstack[ist]; TPZGeoElSide gneigh = cneigh.Reference(); - TPZInterpolatedElement *intelneigh = dynamic_cast(cneigh.Element()); + auto *intelneigh = dynamic_cast(cneigh.Element()); if (!intelneigh) DebugStop(); TPZConnect &con_neigh = intelneigh->Connect(cneigh.Side()); int64_t c_neigh_seqnum = con_neigh.SequenceNumber(); - int con_size = con_neigh.NState() * con_neigh.NShape(); + uint con_size = con_neigh.NState() * con_neigh.NShape(); if (con_size != c_blocksize) DebugStop(); for (int ibl = 0; ibl < con_size; ibl++) { sol.at(block.at(c_neigh_seqnum, 0, ibl, 0)) = sol.at(block.at(c_gelSide_seqnum, 0, ibl, 0)); @@ -991,7 +1008,7 @@ void TPZMHMHDivErrorEstimator::VerifySolutionConsistency(TPZCompMesh* cmesh) { #endif // Checks pressure value on these nodes - TPZInterpolatedElement *intel = dynamic_cast(cneighbour.Element()); + auto *intel = dynamic_cast(cneighbour.Element()); if (!intel) DebugStop(); } } @@ -1082,7 +1099,7 @@ void TPZMHMHDivErrorEstimator::CreateFluxSkeletonElements(TPZCompMesh *flux_mesh if (neighGel->Dimension() != dim) continue; TPZCompEl *neighCel = geltocel[neighGel]; if (!neighCel) DebugStop(); - TPZInterpolatedElement *neighIntel = dynamic_cast(neighCel); + auto *neighIntel = dynamic_cast(neighCel); if (!neighIntel) DebugStop(); // TODO sometimes the order here is 1, I think it should always be 3 for the case I'm running. Need to @@ -1097,7 +1114,7 @@ void TPZMHMHDivErrorEstimator::CreateFluxSkeletonElements(TPZCompMesh *flux_mesh TPZGeoElBC gbc(neighbour, fHDivWrapMatId); int64_t index; TPZCompEl * wrap = flux_mesh->ApproxSpace().CreateCompEl(gbc.CreatedElement(), *flux_mesh, index); - TPZInterpolatedElement *wrapintel = dynamic_cast(wrap); + auto *wrapintel = dynamic_cast(wrap); neighIntel->SetSideOrient(neighbour.Side(), 1); wrapintel->SetSideOrient(gbc.CreatedElement()->NSides()-1,1); int64_t newc_index = wrap->ConnectIndex(0); @@ -1121,7 +1138,7 @@ void TPZMHMHDivErrorEstimator::CreateMultiphysicsInterfaces() { fMultiPhysicsInterfaceMatId = FindFreeMatId(this->GMesh()); } - TPZLagrangeMultiplierCS<> *interfaceMat = new TPZLagrangeMultiplierCS<>(fMultiPhysicsInterfaceMatId, dim - 1, 1); + auto *interfaceMat = new TPZLagrangeMultiplier(fMultiPhysicsInterfaceMatId, dim - 1, 1); fPostProcMesh.InsertMaterialObject(interfaceMat); std::cout << "Created interface material of index " << fMultiPhysicsInterfaceMatId << '\n'; diff --git a/ErrorEstimation/TPZMHMHDivErrorEstimator.h b/ErrorEstimation/TPZMHMHDivErrorEstimator.h index cb286316..af2cc5eb 100644 --- a/ErrorEstimation/TPZMHMHDivErrorEstimator.h +++ b/ErrorEstimation/TPZMHMHDivErrorEstimator.h @@ -5,8 +5,8 @@ // Created by Philippe Devloo on 10/06/18. // -#ifndef TPZMHMHDivErrorEstimator_hpp -#define TPZMHMHDivErrorEstimator_hpp +#ifndef TPZMHMHDIVERRORESTIMATOR_H +#define TPZMHMHDIVERRORESTIMATOR_H #include "TPZHDivErrorEstimator.h" #include "TPZMHMixedMeshControl.h" @@ -85,4 +85,4 @@ class TPZMHMHDivErrorEstimator : public TPZHDivErrorEstimator { void ComputeConnectsNextToSkeleton(std::set& connectList); }; -#endif /* TPZHybridHDivErrorEstimator_hpp */ +#endif \ No newline at end of file diff --git a/Projects/CMakeLists.txt b/Projects/CMakeLists.txt index d1df87cb..b7ea9ac5 100644 --- a/Projects/CMakeLists.txt +++ b/Projects/CMakeLists.txt @@ -4,7 +4,7 @@ add_subdirectory(Tools) # Created executable targets #add_subdirectory(ErrorEstimationH1) add_subdirectory(ErrorEstimationHDiv) -add_subdirectory(HybridH1) +#add_subdirectory(HybridH1) add_subdirectory(ErrorEstimationMHM) # This project is marked as optional since it requires OpenCV external library @@ -16,12 +16,8 @@ endif() # This project is marked as optional since it requires libInterpolate external library option(BUILD_UNISIM_PROJECT "Whether to build estimation on UNISIM mesh" OFF) if(BUILD_UNISIM_PROJECT) - add_subdirectory(UNISIM) + #add_subdirectory(UNISIM) endif() -# This project is marked as optional since it requires libInterpolate external library -option(BUILD_SPE10_PROJECT "Whether to build estimation on SPE10 mesh" OFF) -if(BUILD_SPE10_PROJECT) - add_subdirectory(SPE10) -endif() +add_subdirectory(SPE10) diff --git a/Projects/ErrorEstimationH1/TPZPostProcessError.cpp b/Projects/ErrorEstimationH1/TPZPostProcessError.cpp index fa272504..0f120247 100644 --- a/Projects/ErrorEstimationH1/TPZPostProcessError.cpp +++ b/Projects/ErrorEstimationH1/TPZPostProcessError.cpp @@ -267,8 +267,9 @@ void TPZPostProcessError::ComputeHDivSolution() TPZCompMesh *meshmixed = fMeshVector[1]; int64_t nval = fMeshVector[4]->Solution().Rows(); + TPZFMatrix &sol = fMeshVector[4]->Solution(); for (int64_t i=0; iSolution()(i,0) = 1.; + sol(i,0) = 1.; } int ModelDimension = meshmixed->Dimension(); @@ -371,11 +372,13 @@ void TPZPostProcessError::ComputeElementErrors(TPZVec &elementerrors) } } int64_t extseqcount = numintconnects; + TPZFMatrix &sol = meshweight->Solution(); for (int64_t patch=0; patchConnectVec()[partitionindex].SequenceNumber(); - meshweight->Block()(seqnum,0,0,0) = 1.; + std::pair ind = meshweight->Block().at(seqnum, 0, 0, 0); + sol.at(ind) = 1; } { int64_t nel = fVecVecPatches[color][patch].fElIndices.size(); @@ -491,23 +494,23 @@ void TPZPostProcessError::ComputeElementErrors(TPZVec &elementerrors) direct = 0; an.Assemble(); - - TPZAutoPointer > globmat = an.Solver().Matrix(); + + TPZMatrix &globmat = an.Mesh()->Solution(); for (int64_t p=0; pGetVal(firstlagrangeequation, firstlagrangeequation); + STATE diag = globmat.GetVal(firstlagrangeequation, firstlagrangeequation); diag += 1.; - globmat->Put(firstlagrangeequation, firstlagrangeequation, diag); + globmat.Put(firstlagrangeequation, firstlagrangeequation, diag); } } an.Solve(); - TPZStepSolver &step = dynamic_cast &>(an.Solver()); - std::list singular = step.Singular(); + TPZStepSolver *step = dynamic_cast*>(an.Solver()); + std::list singular = step->Singular(); if (singular.size()) { @@ -695,6 +698,7 @@ int64_t TPZPatch::FirstLagrangeEquation(TPZCompMesh *cmesh) const // Sum the solution stored in fSolution of the second mesh to the fSolution vector void TPZPostProcessError::TransferAndSumSolution(TPZCompMesh *cmesh) { + TPZFMatrix &sol = cmesh->Solution(); int64_t nconnect = cmesh->NConnects(); for (int64_t ic=0; icConnectVec()[ic]; @@ -714,7 +718,7 @@ void TPZPostProcessError::TransferAndSumSolution(TPZCompMesh *cmesh) for (int eq=0; eqSolution()(pos+eq,0); + fSolution(targetpos+eq,0) += sol(pos+eq,0); } } cmesh->Solution().Zero(); diff --git a/Projects/ErrorEstimationHDiv/CMakeLists.txt b/Projects/ErrorEstimationHDiv/CMakeLists.txt index 912b4dec..9805ef89 100644 --- a/Projects/ErrorEstimationHDiv/CMakeLists.txt +++ b/Projects/ErrorEstimationHDiv/CMakeLists.txt @@ -1,15 +1,15 @@ # Project to apply error estimation for H(div) simulations -add_executable(HeteroPermeability main_HeteroPermeability.cpp) -target_link_libraries(HeteroPermeability Tools ErrorEstimationLib) +#add_executable(HeteroPermeability main_HeteroPermeability.cpp) +#target_link_libraries(HeteroPermeability Tools ErrorEstimationLib) add_executable(HybridHDiv main.cpp) target_link_libraries(HybridHDiv Tools ErrorEstimationLib) -add_executable(AdaptivityTest mainAdaptivity.cpp) -target_link_libraries(AdaptivityTest Tools ErrorEstimationLib) +add_executable(AdaptivityHDiv mainAdaptivity.cpp) +target_link_libraries(AdaptivityHDiv Tools ErrorEstimationLib) -add_executable(SteklovProblem mainSteklov.cpp) -target_link_libraries(SteklovProblem Tools ErrorEstimationLib) +#add_executable(SteklovProblem mainSteklov.cpp) +#target_link_libraries(SteklovProblem Tools ErrorEstimationLib) configure_file(../Meshes/Quad.msh Quad.msh COPYONLY) configure_file(../Meshes/Cube.msh Cube.msh COPYONLY) diff --git a/Projects/ErrorEstimationHDiv/main.cpp b/Projects/ErrorEstimationHDiv/main.cpp index 4a31cb56..c74a91f9 100644 --- a/Projects/ErrorEstimationHDiv/main.cpp +++ b/Projects/ErrorEstimationHDiv/main.cpp @@ -7,12 +7,12 @@ #include "ProblemConfig.h" #include "TPZBFileStream.h" -#include "TPZGmshReader.h" #include "TPZHDivErrorEstimator.h" #include "TPZHybridizeHDiv.h" #include "TPZMultiphysicsCompMesh.h" #include "TPZRefPatternDataBase.h" #include "Tools.h" +#include "TPZBndCond.h" #include "pzlog.h" #include "tpzgeoelrefpattern.h" #include @@ -54,7 +54,7 @@ void RunSingularProblemHDiv() { config.bcmaterialids.insert(-1); config.makepressurecontinuous = true; - int nRef = 2; + int nRef = 3; config.ndivisions = 0; config.gmesh = CreateLShapeGeoMesh(nRef); @@ -65,7 +65,7 @@ void RunSingularProblemHDiv() { TPZHybridizeHDiv hybridizer; TPZMultiphysicsCompMesh *hybridMesh = CreateHybridCompMesh(config, hybridizer); - Tools::SolveHybridProblem(hybridMesh, hybridizer.fInterfaceMatid, config, false); + Tools::SolveHybridProblem(hybridMesh, hybridizer.fInterfaceMatid, config, true); // Reconstruct potential and estimate error EstimateError(config, hybridMesh, hybridizer); @@ -77,9 +77,9 @@ void RunSingularProblemHDiv() { TPZMultiphysicsCompMesh *CreateHybridCompMesh(const ProblemConfig &config, TPZHybridizeHDiv &hybridizer) { - TPZMultiphysicsCompMesh *cmesh_HDiv = Tools::CreateHDivMesh(config); // Hdiv x L2 + TPZMultiphysicsCompMesh *cmesh_HDiv = Tools::CreateMixedMesh(config); // Hdiv x L2 -#ifdef ERRORESTIMATION_DEBUG +#ifdef PZDEBUG { std::ofstream out("MixedMesh.txt"); cmesh_HDiv->Print(out); diff --git a/Projects/ErrorEstimationHDiv/mainAdaptivity.cpp b/Projects/ErrorEstimationHDiv/mainAdaptivity.cpp index c8120a29..d7e91ef2 100644 --- a/Projects/ErrorEstimationHDiv/mainAdaptivity.cpp +++ b/Projects/ErrorEstimationHDiv/mainAdaptivity.cpp @@ -2,34 +2,23 @@ // Created by Gustavo Batistela on 12/07/19. // -#include "TPZRefPatternDataBase.h" -#include "TPZGmshReader.h" - -#include "tpzarc3d.h" -#include "tpzgeoblend.h" -#include "TPZGeoLinear.h" - -#include -#include - #include "ProblemConfig.h" #include "TPZHDivErrorEstimator.h" +#include "TPZRefPatternDataBase.h" #include "Tools.h" -//#include "pzelchdiv.h" - -bool readGeoMeshFromFile = false; -bool postProcessWithHDiv = false; -int refinementSteps = 1; +void RunAdaptivitySuite(int refinementSteps); +void RunUniformRefinementSuite(int refinementSteps); -TPZGeoMesh *CreateLShapeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &mhmIndexes); -void TracingTriangleBug(TPZMultiphysicsCompMesh* multiphysics); +int main() { -int main(int argc, char* argv[]) { + constexpr int adaptivitySteps = 8; + constexpr int uniformSteps = 6; + RunAdaptivitySuite(adaptivitySteps); + RunUniformRefinementSuite(uniformSteps); +} -#ifdef LOG4CXX - InitializePZLOG(); -#endif +void RunAdaptivitySuite(int refinementSteps) { // Initializing uniform refinements for reference elements gRefDBase.InitializeUniformRefPattern(EOned); @@ -37,182 +26,129 @@ int main(int argc, char* argv[]) { gRefDBase.InitializeUniformRefPattern(ETriangle); // Creates geometric mesh - - TPZGeoMesh *gmeshOriginal = nullptr; // CreateLShapeMesh(bcIDs);//CreateLShapeMesh(bcIDs); + TPZGeoMesh *gmeshOriginal; ProblemConfig config; config.porder = 1; - config.hdivmais = 1; + config.hdivmais = 2; config.dimension = 2; config.makepressurecontinuous = true; config.exact = new TLaplaceExample1; - config.exact.operator*().fExact = TLaplaceExample1::ESinMark; + config.exact.operator*().fExact = TLaplaceExample1::EBoundaryLayer; - config.dir_name = "AdaptivityLShape"; - config.problemname = "ESinSinMark"; + config.dir_name = "HDivAdaptivityEquations"; + config.problemname = "Adaptivity"; + { + std::string command = "mkdir -p " + config.dir_name; + system(command.c_str()); + } - std::string command = "mkdir -p " + config.dir_name; - system(command.c_str()); + TPZManVector bcids(4, -1); + gmeshOriginal = Tools::CreateGeoMesh(3, bcids); + config.materialids.insert(1); + config.bcmaterialids.insert(-1); - if (readGeoMeshFromFile) { - std::string meshfilename = "../LMesh.msh"; //"../LMesh3.msh"; + Tools::DivideLowerDimensionalElements(gmeshOriginal); - TPZGmshReader gmsh; - gmsh.GetDimNamePhysical()[1]["dirichlet"] = 2; - gmsh.GetDimNamePhysical()[2]["domain"] = 1; - gmeshOriginal = gmsh.GeometricGmshMesh(meshfilename); - gmsh.PrintPartitionSummary(std::cout); - config.materialids.insert(1); - config.bcmaterialids.insert(2); + for (int iStep = 0; iStep < refinementSteps; iStep++) { - } -else + config.gmesh = new TPZGeoMesh(*gmeshOriginal); + config.adaptivityStep = iStep; - { - //TPZManVector bcIDs(8, -1); - TPZManVector bcids(8,-1); - bcids[1] = -1; - //constants for Robin boundary conditions - // sigma.n=Km(u-u_d)-g - //Particular cases: 1) Km=0---> Neumann, 2) Km=infinity-->Dirichlet - //config.coefG = 0.;//nao passar mais isso - config.Km = 1.e12;//pow(10,2); - - gmeshOriginal = Tools::CreateLShapeMesh(bcids);//CreateGeoMesh(1, bcIDs);// - - - config.materialids.insert(1); - config.bcmaterialids.insert(-1); // dirichlet - //config.bcmaterialids.insert(-2); // neumann - //config.bcmaterialids.insert(-3); // Robin + TPZMultiphysicsCompMesh *mixedCompMesh = Tools::CreateMixedMesh(config); // Hdiv x L2 + mixedCompMesh->InitializeBlock(); + Tools::SolveMixedProblem(mixedCompMesh, config); + config.cmesh = mixedCompMesh; - - } - - gmeshOriginal->SetDimension(config.dimension); - config.gmesh = gmeshOriginal; - Tools::UniformRefinement(2, 2 , config.gmesh) ; - Tools::DivideLowerDimensionalElements(config.gmesh); - - - - TPZManVector meshvec_HDiv(2, 0); - - TPZMultiphysicsCompMesh* cmesh_HDiv = Tools::CreateHDivMesh(config); //Hdiv x L2 - cmesh_HDiv->InitializeBlock(); - Tools::SolveMixedProblem(cmesh_HDiv, config); - - meshvec_HDiv = cmesh_HDiv->MeshVector(); - - - for (int iSteps = 0; iSteps < refinementSteps; iSteps++) { - - - config.adaptivityStep = iSteps; - - // UniformRefinement(iSteps, gmeshOriginal); - - #ifdef ERRORESTIMATION_DEBUG - { - std::ofstream out("gmeshToSolve.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(gmeshOriginal, out); - } - #endif - - - - - TPZManVector meshvec_HDiv(2, 0); - - TPZMultiphysicsCompMesh* cmesh_HDiv = nullptr; - - - cmesh_HDiv = Tools::CreateHDivMesh(config);//Hdiv x L2 - cmesh_HDiv->InitializeBlock(); - #ifdef ERRORESTIMATION_DEBUG2 - { - std::ofstream out("MultiPhysicsMesh.txt"); - cmesh_HDiv->Print(out); - std::ofstream outvtk("MultiPhysicsMesh.vtk"); - - TPZVTKGeoMesh::PrintCMeshVTK(cmesh_HDiv,outvtk); - - - } - #endif - - meshvec_HDiv = cmesh_HDiv->MeshVector(); - - //cria malha hibrida - TPZHybridizeHDiv hybrid; - auto HybridMesh = hybrid.Hybridize(cmesh_HDiv); - HybridMesh->CleanUpUnconnectedNodes(); - HybridMesh->AdjustBoundaryElements(); - - delete cmesh_HDiv; - delete meshvec_HDiv[0]; - delete meshvec_HDiv[1]; - - cmesh_HDiv = (HybridMesh);//malha hribrida - meshvec_HDiv[0] = (HybridMesh)->MeshVector()[0];//malha Hdiv - meshvec_HDiv[1] = (HybridMesh)->MeshVector()[1];//malha L2 - - Tools::SolveHybridProblem(cmesh_HDiv, hybrid.fInterfaceMatid, config,false); - - - - //reconstroi potencial e calcula o erro + // Estimate error and run adaptive process { bool postProcWithHDiv = false; - TPZHDivErrorEstimator HDivEstimate(*cmesh_HDiv, postProcWithHDiv); + TPZHDivErrorEstimator HDivEstimate(*mixedCompMesh, postProcWithHDiv); HDivEstimate.SetAnalyticSolution(config.exact); - + HDivEstimate.SetAdaptivityStep(iStep); HDivEstimate.PotentialReconstruction(); TPZManVector elementerrors; TPZManVector errorvec; - std::string vtkPath = "adaptivity_error_results.vtk"; - HDivEstimate.ComputeErrors(errorvec, elementerrors, vtkPath); + std::stringstream outVTK; + outVTK << config.dir_name << "/" << config.problemname << "-Errors.vtk"; + auto stringVTK = outVTK.str(); + HDivEstimate.ComputeErrors(errorvec, elementerrors, stringVTK); Tools::hAdaptivity(HDivEstimate.PostProcMesh(), gmeshOriginal, config); - #ifdef ERRORESTIMATION_DEBUG - { - std::ofstream out("gmeshAdapty.vtk"); - TPZVTKGeoMesh::PrintGMeshVTK(gmeshOriginal, out); - } - #endif - } - - // return 0; - + + std::stringstream outTXT; + outTXT << config.dir_name << "/" << config.problemname << "-Errors-Step" << config.adaptivityStep << ".txt"; + std::ofstream fileTXT(outTXT.str()); + Tools::PrintErrors(fileTXT, config, errorvec); + } } } -TPZGeoMesh *CreateLShapeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &mhmIndexes) { +void RunUniformRefinementSuite(int refinementSteps) { - TPZVec bcIDs(8, -1); - TPZGeoMesh *gmesh = Tools::CreateQuadLShapeMesh(bcIDs); - gmesh->SetDimension(2); - gmesh->BuildConnectivity(); + // Initializing uniform refinements for reference elements + gRefDBase.InitializeUniformRefPattern(EOned); + gRefDBase.InitializeUniformRefPattern(EQuadrilateral); + gRefDBase.InitializeUniformRefPattern(ETriangle); - Tools::UniformRefinement(nCoarseRef, gmesh); - Tools::DivideLowerDimensionalElements(gmesh); + ProblemConfig config; - int64_t nElem = gmesh->NElements(); - for (int64_t i = 0; i < nElem; i++) { - TPZGeoEl *gel = gmesh->Element(i); - if (gel->Dimension() != gmesh->Dimension() || gel->NSubElements() > 0) continue; - mhmIndexes.Push(i); - } + config.porder = 1; + config.hdivmais = 2; - Tools::UniformRefinement(nInternalRef, gmesh); - Tools::DivideLowerDimensionalElements(gmesh); + config.dimension = 2; + config.makepressurecontinuous = true; - for (int64_t i = 0; i < mhmIndexes.size(); i++) { - std::cout << mhmIndexes[i] << '\n'; + config.exact = new TLaplaceExample1; + config.exact.operator*().fExact = TLaplaceExample1::EBoundaryLayer; + + config.dir_name = "HDivAdaptivityEquations"; + config.problemname = "UnifRef"; + { + std::string command = "mkdir -p " + config.dir_name; + system(command.c_str()); } - std::cout << '\n'; - return gmesh; + for (int iStep = 0; iStep < refinementSteps; iStep++) { + + // Creates geometric mesh + TPZGeoMesh *gmeshOriginal; + + TPZManVector bcids(4, -1); + int nelem = 3 * static_cast(pow(2.0, iStep)); + gmeshOriginal = Tools::CreateGeoMesh(nelem, bcids); + config.materialids.insert(1); + config.bcmaterialids.insert(-1); + + Tools::DivideLowerDimensionalElements(gmeshOriginal); + + config.gmesh = new TPZGeoMesh(*gmeshOriginal); + config.adaptivityStep = iStep; + + TPZMultiphysicsCompMesh *mixedCompMesh = Tools::CreateMixedMesh(config); // Hdiv x L2 + mixedCompMesh->InitializeBlock(); + Tools::SolveMixedProblem(mixedCompMesh, config); + config.cmesh = mixedCompMesh; + + { + bool postProcessWithHDiv = false; + TPZHDivErrorEstimator HDivEstimate(*mixedCompMesh, postProcessWithHDiv); + HDivEstimate.SetAnalyticSolution(config.exact); + HDivEstimate.SetAdaptivityStep(iStep); + HDivEstimate.PotentialReconstruction(); + TPZManVector elementerrors; + TPZManVector errorvec; + std::stringstream outVTK; + outVTK << config.dir_name << "/" << config.problemname << "-Errors.vtk"; + auto stringVTK = outVTK.str(); + HDivEstimate.ComputeErrors(errorvec, elementerrors, stringVTK); + + std::stringstream outTXT; + outTXT << config.dir_name << "/" << config.problemname << "-Errors-Step" << config.adaptivityStep << ".txt"; + std::ofstream fileTXT(outTXT.str()); + Tools::PrintErrors(fileTXT, config, errorvec); + } + } } diff --git a/Projects/ErrorEstimationMHM/CMakeLists.txt b/Projects/ErrorEstimationMHM/CMakeLists.txt index c8db4df8..c01e420b 100644 --- a/Projects/ErrorEstimationMHM/CMakeLists.txt +++ b/Projects/ErrorEstimationMHM/CMakeLists.txt @@ -4,3 +4,6 @@ target_link_libraries(MHM_Estimation Tools ErrorEstimationLib) add_executable(MHM_Journal main_Journal.cpp) target_link_libraries(MHM_Journal Tools ErrorEstimationLib) + +add_executable(MHM_RefSol main_RefSolValidation.cpp) +target_link_libraries(MHM_RefSol Tools ErrorEstimationLib) diff --git a/Projects/ErrorEstimationMHM/main_Journal.cpp b/Projects/ErrorEstimationMHM/main_Journal.cpp index 21ae6494..f429106d 100644 --- a/Projects/ErrorEstimationMHM/main_Journal.cpp +++ b/Projects/ErrorEstimationMHM/main_Journal.cpp @@ -1,68 +1,149 @@ // // Created by Gustavo Batistela on 3/31/21. // - +#include "ProblemConfig.h" +#include "TPZMHMHDivErrorEstimator.h" +#include "ToolsMHM.h" #include "pzgmesh.h" +#include #include
 #include 
+#include 
+#include 
 #include 
-//#include 
-#include 
+#include 
+#include 
+#include 
 #include 
-
-void RunSmoothProblem(int nCoarseDiv, int nInternalRef);
+#include 
+#include 
+
+struct ErrorResult {
+    ErrorResult(const std::string &s, const int nIntRef, const int nCoarseDiv, const int kOrder, const int nOrder,
+                const REAL estimatedError, const REAL exactError, const REAL ieff) {
+        problem_name = s;
+        n_internal_ref = nIntRef;
+        n_coarse_div = nCoarseDiv;
+        k_order = kOrder;
+        n_order = nOrder;
+        estimated_error = estimatedError;
+        exact_error = exactError;
+        effectivity_index = ieff;
+    }
+    std::string problem_name = "NULL";
+    int n_internal_ref = -1;
+    int n_coarse_div = -1;
+    int k_order = -1;
+    int n_order = -1;
+    REAL estimated_error = -1.;
+    REAL exact_error = -1.;
+    REAL effectivity_index = -1.;
+};
+
+std::vector result_vec;
+
+const int porder = 1;
+const int hdivmais = 4;
+std::string dir_name = "Journal_k" + std::to_string(porder) + "n" + std::to_string(hdivmais);
+
+void RunSmoothProblem(int nCoarseDiv, int nInternalRef, int k, int n);
 void RunNonConvexProblem();
 void RunHighGradientProblem(int nCoarseDiv, int nInternalRef);
-void RunInnerSingularityProblem(int nCoarseDiv, int nInternalRef);
+void RunInnerSingularityProblem(int nCoarseDiv, int nInternalRef, int k, int n);
 void RunPeriodicPermProblem(int nCoarseDiv, int nInternalRef);
+void RunHighGradientAdaptivityProblem(int n_divisions);
 
 TPZGeoMesh *CreateQuadGeoMesh(int nCoarseDiv, int nInternalRef);
 TPZGeoMesh *CreateLShapeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &mhmIndexes);
 
 void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfig &config);
-void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, int nInternalRef,
-                       bool definePartitionByCoarseIndex, TPZManVector& mhmIndexes);
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes);
+
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes, const std::vector &refLevelPerSubdomain);
+
 void CreateMHMCompMeshHeteroPerm(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, int nInternalRef,
                                  bool definePartitionByCoarseIndex, TPZManVector& mhmIndexes);
 
 void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config);
 
 void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm);
+void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm, TPZMHMHDivErrorEstimator &estimator);
 
-void MHMAdaptivity(TPZMHMixedMeshControl *mhm, TPZGeoMesh* gmeshToRefine, ProblemConfig& config);
+void MHMAdaptivity(TPZMHMixedMeshControl *mhm, ProblemConfig &config, TPZCompMesh *postProcMesh,
+                   std::vector &refLevelPerSubdomain);
 
 void CreateMHMCompMeshPermFunction(TPZMHMixedMeshControl &mhm);
-void PeriodicPermeabilityFunction(const TPZVec &coord, TPZVec &res, TPZFMatrix &res_mat);
 void PeriodicProblemForcingFunction(const TPZVec  &pt, TPZVec  &result);
 
+void PrintLatexGraphs(std::ostream & out);
+
+void RunSmoothProblemSuite(const std::set &nCoarseDiv, const std::set &nInternalRef,
+                           const std::set &kOrder, const std::set &nOrder);
+void RunInnerSingularityProblemSuite(const std::set &nCoarseDiv, const std::set &nInternalRef,
+                                     const std::set &kOrder, const std::set &nOrder);
+
+
+
 int main() {
+#ifdef PZ_LOG
     TPZLogger::InitializePZLOG();
+#endif
     gRefDBase.InitializeAllUniformRefPatterns();
 
-    const std::set nCoarseDiv = {3, 4, 5, 6};
-    const std::set nInternalRef = {0, 1, 2, 3};
-    for (const auto coarse_div : nCoarseDiv) {
-        for (const auto internal_ref : nInternalRef) {
-            RunSmoothProblem(coarse_div, internal_ref);
-            //RunInnerSingularityProblem(coarse_div, internal_ref);
-            //RunHighGradientProblem(coarse_div, internal_ref);
-            //RunPeriodicPermProblem(coarse_div, internal_ref);
-        }
-    }
-    RunNonConvexProblem();
+    const std::set nCoarseDiv = {8};
+    const std::set nInternalRef = {3};
+    const std::set kOrder = {1};
+    const std::set nOrder = {2};
+
+    RunSmoothProblemSuite(nCoarseDiv, nInternalRef, kOrder, nOrder);
+
+    //RunInnerSingularityProblemSuite(nCoarseDiv, nInternalRef, kOrder, nOrder);
+
+    //RunHighGradientAdaptivityProblem(4);
+
+    //RunNonConvexProblem();
+    //std::ofstream out_latex("LatexGraphs.txt", std::ios::trunc);
+    //PrintLatexGraphs(out_latex);
 
     return 0;
 }
 
-void RunSmoothProblem(const int nCoarseDiv, const int nInternalRef) {
+void RunInnerSingularityProblemSuite(const std::set &nCoarseDiv, const std::set &nInternalRef,
+                                     const std::set &kOrder, const std::set &nOrder) {
+    for (const auto k : kOrder) {
+        for (const auto n : nOrder) {
+            for (const auto internal_ref : nInternalRef) {
+                for (const auto coarse_div : nCoarseDiv) {
+                    RunInnerSingularityProblem(coarse_div, internal_ref, k, n);
+                }
+            }
+        }
+    }
+}
+void RunSmoothProblemSuite(const std::set &nCoarseDiv, const std::set &nInternalRef,
+                           const std::set &kOrder, const std::set &nOrder) {
+    for (const auto k : kOrder) {
+        for (const auto n : nOrder) {
+            for (const auto internal_ref : nInternalRef) {
+                for (const auto coarse_div : nCoarseDiv) {
+                    RunSmoothProblem(coarse_div, internal_ref, k, n);
+                }
+            }
+        }
+    }
+}
+
+void RunSmoothProblem(const int nCoarseDiv, const int nInternalRef, const int k, const int n) {
     ProblemConfig config;
     config.dimension = 2;
     config.exact = new TLaplaceExample1;
     config.exact.operator*().fExact = TLaplaceExample1::ESinSin;
     config.problemname = "Smooth";
-    config.dir_name = "Journal";
-    config.porder = 1;
-    config.hdivmais = 2;
+    config.dir_name = dir_name;
+    config.porder = k;
+    config.hdivmais = n;
     config.materialids.insert(1);
     config.bcmaterialids.insert(-1);
     config.makepressurecontinuous = true;
@@ -84,9 +165,10 @@ void RunSmoothProblem(const int nCoarseDiv, const int nInternalRef) {
     TPZManVector coarseIndexes;
     ComputeCoarseIndices(config.gmesh, coarseIndexes);
     bool definePartitionByCoarseIndexes = true;
-    CreateMHMCompMesh(mhm, config, nInternalRef, definePartitionByCoarseIndexes, coarseIndexes);
+    CreateMHMCompMesh(mhm, config, definePartitionByCoarseIndexes, coarseIndexes);
 
     SolveMHMProblem(mhm, config);
+    TPZCompMesh *cmesh = nullptr;
     EstimateError(config, mhm);
 }
 
@@ -96,9 +178,9 @@ void RunHighGradientProblem(const int nCoarseDiv, const int nInternalRef) {
     config.exact = new TLaplaceExample1;
     config.exact.operator*().fExact = TLaplaceExample1::EBoundaryLayer;
     config.problemname = "HighGradient";
-    config.dir_name = "Journal";
-    config.porder = 1;
-    config.hdivmais = 2;
+    config.dir_name = dir_name;
+    config.porder = porder;
+    config.hdivmais = hdivmais;
     config.materialids.insert(1);
     config.bcmaterialids.insert(-1);
     config.makepressurecontinuous = true;
@@ -120,7 +202,7 @@ void RunHighGradientProblem(const int nCoarseDiv, const int nInternalRef) {
     TPZManVector coarseIndexes;
     ComputeCoarseIndices(config.gmesh, coarseIndexes);
     bool definePartitionByCoarseIndexes = true;
-    CreateMHMCompMesh(mhm, config, nInternalRef, definePartitionByCoarseIndexes, coarseIndexes);
+    CreateMHMCompMesh(mhm, config, definePartitionByCoarseIndexes, coarseIndexes);
 
     SolveMHMProblem(mhm, config);
     EstimateError(config, mhm);
@@ -132,9 +214,9 @@ void RunNonConvexProblem() {
     config.exact = new TLaplaceExample1;
     config.exact.operator*().fExact = TLaplaceExample1::ESinSin;
     config.problemname = "NonConvex";
-    config.dir_name = "Journal";
-    config.porder = 1;
-    config.hdivmais = 3;
+    config.dir_name = dir_name;
+    config.porder = porder;
+    config.hdivmais = hdivmais;
     config.materialids.insert(1);
     config.bcmaterialids.insert(-1);
     config.makepressurecontinuous = true;
@@ -157,21 +239,22 @@ void RunNonConvexProblem() {
 
     auto *mhm = new TPZMHMixedMeshControl(config.gmesh);
     bool definePartitionByCoarseIndexes = false;
-    CreateMHMCompMesh(mhm, config, nInternalRef, definePartitionByCoarseIndexes, coarseIndexes);
+    CreateMHMCompMesh(mhm, config, definePartitionByCoarseIndexes, coarseIndexes);
 
     SolveMHMProblem(mhm, config);
+    TPZCompMesh *cmesh = nullptr;
     EstimateError(config, mhm);
 }
 
-void RunInnerSingularityProblem(const int nCoarseDiv, const int nInternalRef) {
+void RunInnerSingularityProblem(const int nCoarseDiv, const int nInternalRef, int k, int n) {
     ProblemConfig config;
     config.dimension = 2;
     config.exact = new TLaplaceExample1;
     config.exact.operator*().fExact = TLaplaceExample1::ESteklovNonConst;
     config.problemname = "InnerSingularity";
-    config.dir_name = "Journal";
-    config.porder = 1;
-    config.hdivmais = 2;
+    config.dir_name = dir_name;
+    config.porder = k;
+    config.hdivmais = n;
     config.materialids = {1, 2};
     config.bcmaterialids.insert(-1);
     config.makepressurecontinuous = true;
@@ -228,7 +311,9 @@ void RunInnerSingularityProblem(const int nCoarseDiv, const int nInternalRef) {
     CreateMHMCompMeshHeteroPerm(mhm, config, nInternalRef, definePartitionByCoarseIndexes, coarseIndexes);
 
     SolveMHMProblem(mhm, config);
+    TPZCompMesh *cmesh = nullptr;
     EstimateError(config, mhm);
+
 }
 
 void RunPeriodicPermProblem(const int nCoarseDiv, const int nInternalRef) {
@@ -238,9 +323,9 @@ void RunPeriodicPermProblem(const int nCoarseDiv, const int nInternalRef) {
     //config.exact = new TLaplaceExample1;
     //config.exact.operator*().fExact = TLaplaceExample1::EBoundaryLayer;
     config.problemname = "PeriodicPerm";
-    config.dir_name = "Journal";
-    config.porder = 1;
-    config.hdivmais = 2;
+    config.dir_name = dir_name;
+    config.porder = porder;
+    config.hdivmais = hdivmais;
     config.materialids.insert(1);
     config.bcmaterialids.insert(-1);
     config.makepressurecontinuous = true;
@@ -389,7 +474,45 @@ TPZGeoMesh *CreateLShapeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &mhmIndexes, const std::vector &refLevelPerSubdomain) {
+
+    if (definePartitionByCoarseIndex) {
+        mhm->DefinePartitionbyCoarseIndices(mhmIndexes);
+    } else {
+        mhm->DefinePartition(mhmIndexes);
+    }
+
+    // Indicate material indices to the MHM control structure
+    mhm->fMaterialIds = config.materialids;
+    mhm->fMaterialBCIds = config.bcmaterialids;
+
+    // Insert the material objects in the multiphysics mesh
+    InsertMaterialsInMHMMesh(*mhm, config);
+
+    // General approximation order settings
+    mhm->SetInternalPOrder(config.porder);
+    mhm->SetSkeletonPOrder(config.porder);
+    mhm->SetHdivmaismaisPOrder(config.hdivmais);
+
+    const auto interfaces = mhm->GetInterfaces();
+    if (!refLevelPerSubdomain.empty()) {
+        for (auto interface : interfaces) {
+            if (interface.first == interface.second.first || interface.first == interface.second.second) continue;
+            auto right_ref_level = refLevelPerSubdomain[interface.second.first];
+            auto left_ref_level = refLevelPerSubdomain[interface.second.second];
+            const auto n_divisions = std::max(right_ref_level, left_ref_level);
+            mhm->DivideSkeletonElement(interface.first, n_divisions);
+        }
+    }
+
+    mhm->DivideBoundarySkeletonElements();
+    // Creates MHM mesh
+    bool substructure = true;
+    mhm->BuildComputationalMesh(substructure);
+}
+
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config,
                        bool definePartitionByCoarseIndex, TPZManVector& mhmIndexes) {
 
     if (definePartitionByCoarseIndex) {
@@ -423,13 +546,13 @@ void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config) {
     TPZAutoPointer cmesh = mhm->CMesh();
 
     bool shouldrenumber = true;
-    TPZAnalysis an(cmesh, shouldrenumber);
+    TPZLinearAnalysis an(cmesh, shouldrenumber);
 
 #ifdef PZ_USING_MKL
-    TPZSymetricSpStructMatrix strmat(cmesh.operator->());
-    strmat.SetNumThreads(8);
+    TPZSSpStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(4);
 #else
-    TPZSkylineStructMatrix strmat(cmesh.operator->());
+    TPZSkylineStructMatrix strmat(cmesh.operator->());
     strmat.SetNumThreads(0);
 #endif
 
@@ -445,7 +568,6 @@ void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config) {
     std::cout << "Finished\n";
     an.LoadSolution(); // compute internal dofs
 
-    // TODO: Phil, this is the part that needs a fix!
     TPZMFSolutionTransfer transfer;
     transfer.BuildTransferData(cmesh.operator->());
     transfer.TransferFromMultiphysics();
@@ -472,7 +594,7 @@ void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config) {
     plotname << config.dir_name << "/" << config.problemname << "-" << config.ndivisions << "-" << config.ninternalref
              << "-Results.vtk";
     an.DefineGraphMesh(cmesh->Dimension(), scalnames, vecnames, plotname.str());
-    an.PostProcess(resolution, cmesh->Dimension());
+    //an.PostProcess(resolution, cmesh->Dimension());
 }
 
 void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm) {
@@ -504,6 +626,16 @@ void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm) {
         std::ofstream file(fileName, std::ios::app);
         Tools::PrintErrors(file, config, errors);
     }
+
+    result_vec.emplace_back(config.problemname,
+                            config.ninternalref,
+                            config.ndivisions,
+                            config.porder,
+                            config.hdivmais,
+                            errors[3],
+                            errors[2],
+                            sqrt(errors[4] + errors[3]) / sqrt(errors[2])
+                            );
 }
 
 void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfig &config) {
@@ -512,23 +644,33 @@ void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfi
     int dim = control.GMesh()->Dimension();
     cmesh.SetDimModel(dim);
 
-    TPZMixedPoisson *mat = new TPZMixedPoisson(1, dim);
+    auto *mat = new TPZMixedDarcyFlow(1, dim);
 
-    TPZFMatrix K(3, 3, 0), invK(3, 3, 0);
-    K.Identity();
-    invK.Identity();
+    auto ff_lambda = [config](const TPZVec &loc, TPZVec &result) {
+        config.exact.operator*().ForcingFunction()->Execute(loc, result);
+    };
+    auto exact_sol_lambda = [config](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+        config.exact.operator*().Exact()->Execute(loc, result, deriv);
+    };
 
-    mat->SetExactSol(config.exact.operator*().Exact());
-    mat->SetForcingFunction(config.exact.operator*().ForcingFunction());
-    mat->SetPermeabilityTensor(K, invK);
+    mat->SetForcingFunction(ff_lambda, 5);
+    mat->SetExactSol(exact_sol_lambda, 5);
+    mat->SetConstantPermeability(1);
 
     cmesh.InsertMaterialObject(mat);
 
     for (auto matid : config.bcmaterialids) {
-        TPZFNMatrix<1, REAL> val1(1, 1, 0.), val2(1, 1, 0.);
+        TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+        TPZManVector val2(1, 0.);
         int bctype = 0;
-        TPZBndCond *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-        bc->TPZMaterial::SetForcingFunction(config.exact.operator*().Exact());
+        auto *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+
+        auto ff_bc_lambda = [config](const TPZVec &loc,
+                                     TPZVec &rhsVal,
+                                     TPZFMatrix &matVal) {
+            config.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_bc_lambda);
         cmesh.InsertMaterialObject(bc);
     }
 }
@@ -552,39 +694,40 @@ void CreateMHMCompMeshHeteroPerm(TPZMHMixedMeshControl *mhm, const ProblemConfig
     int dim = mhm->GMesh()->Dimension();
     cmesh.SetDimModel(dim);
 
-    TPZMixedPoisson *mat = new TPZMixedPoisson(1, dim);
+    auto *mat = new TPZMixedDarcyFlow(1, dim);
 
-    TPZFMatrix K(3, 3, 0), invK(3, 3, 0);
-    K.Identity();
-    K *= 5.;
-    invK.Identity();
-    invK *= 1./5.;
+    auto ff_lambda = [config](const TPZVec &loc, TPZVec &result) {
+        config.exact.operator*().ForcingFunction()->Execute(loc, result);
+    };
+    auto exact_sol_lambda = [config](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+        config.exact.operator*().Exact()->Execute(loc, result, deriv);
+    };
 
-    mat->SetExactSol(config.exact.operator*().Exact());
-    mat->SetForcingFunction(config.exact.operator*().ForcingFunction());
-    mat->SetPermeabilityTensor(K, invK);
+    mat->SetForcingFunction(ff_lambda, 5);
+    mat->SetExactSol(exact_sol_lambda, 5);
+    mat->SetConstantPermeability(5);
 
     cmesh.InsertMaterialObject(mat);
 
-    TPZMixedPoisson *mat2 = new TPZMixedPoisson(2, dim);
-
-    TPZFMatrix K2(3, 3, 0), invK2(3, 3, 0);
-    K2.Identity();
-    K2 *= 1.;
-    invK2.Identity();
-    invK2 *= 1./1.;
-
-    mat2->SetExactSol(config.exact.operator*().Exact());
-    mat2->SetForcingFunction(config.exact.operator*().ForcingFunction());
-    mat2->SetPermeabilityTensor(K2, invK2);
+    auto *mat2 = new TPZMixedDarcyFlow(2, dim);
+    mat2->SetForcingFunction(ff_lambda, 5);
+    mat2->SetExactSol(exact_sol_lambda, 5);
+    mat2->SetConstantPermeability(1);
 
     cmesh.InsertMaterialObject(mat2);
 
     for (auto matid : config.bcmaterialids) {
-        TPZFNMatrix<1, REAL> val1(1, 1, 0.), val2(1, 1, 0.);
+        TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+        TPZManVector val2(1, 0.);
         int bctype = 0;
-        TPZBndCond *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-        bc->TPZMaterial::SetForcingFunction(config.exact.operator*().Exact());
+        auto *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+
+        auto ff_bc_lambda = [config](const TPZVec &loc,
+                                     TPZVec &rhsVal,
+                                     TPZFMatrix &matVal) {
+            config.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_bc_lambda);
         cmesh.InsertMaterialObject(bc);
     }
 
@@ -619,12 +762,26 @@ void CreateMHMCompMeshPermFunction(TPZMHMixedMeshControl &mhm) {
 
     // Insert the material objects in the multiphysics mesh
     TPZCompMesh *cmesh = mhm.CMesh().operator->();
-    auto *mix = new TPZMixedPoisson(1, cmesh->Dimension());
-    TPZAutoPointer> perm_function = new TPZDummyFunction(&PeriodicPermeabilityFunction, 3);
+    auto *mix = new TPZMixedDarcyFlow(1, cmesh->Dimension());
+
+    std::function &coord)> perm_function = [](const TPZVec &coord) {
+        constexpr auto epsilon = 0.04;
+        constexpr auto P = 1.8;
+        const auto x = coord[0];
+        const auto y = coord[1];
+
+        REAL term_1 = 2 + P * cos(2 * M_PI * (x - 0.5) / epsilon);
+        REAL term_2 = 2 + P * cos(2 * M_PI * (y - 0.5) / epsilon);
+
+        auto perm = 1 / (term_1 * term_2);
+        return perm;
+    };
+
     mix->SetPermeabilityFunction(perm_function);
     mix->SetForcingFunction(PeriodicProblemForcingFunction, 1);
 
-    TPZFNMatrix<1, REAL> val1(1, 1, 0.), val2(1, 1, 0.);
+    TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+    TPZManVector val2(1, 0.);
     constexpr int dirichlet_bc = 0;
     TPZBndCond *pressure_left = mix->CreateBC(mix, -1, dirichlet_bc, val1, val2);
 
@@ -645,22 +802,287 @@ void CreateMHMCompMeshPermFunction(TPZMHMixedMeshControl &mhm) {
     mhm.BuildComputationalMesh(substructure);
 }
 
-void PeriodicPermeabilityFunction(const TPZVec &coord, TPZVec &res, TPZFMatrix &res_mat) {
+void PeriodicProblemForcingFunction(const TPZVec &pt, TPZVec &result) { result[0] = -1; }
+
+void PrintLatexGraphs(std::ostream & out) {
+
+    std::stringstream latex_text;
+    if (result_vec.empty()) DebugStop();
+
+    auto it = result_vec.begin();
+    while (it != result_vec.end()) {
+        std::string current_problem = it->problem_name;
+        latex_text << "\\section{" << it->problem_name << "}\n";
+        while (it->problem_name == current_problem) {
+
+            int previous_k = it->k_order;
+            while (it->k_order == previous_k) {
+                int previous_n = it->n_order;
+                while (it->n_order == previous_n) {
+                    latex_text << "  \\subsection{$k = " << it->k_order << ", n = " << it->n_order << "$}\n";
+
+                    std::stringstream error_graph, ieff_graph, estimated_error_pts, exact_error_pts,
+                        estimated_error_legend, exact_error_legend, ieff_error_legend;
+
+                    error_graph << "    \\begin{figure}[ht!]\n"
+                                   "      \\centering\n"
+                                   "      \\begin{tikzpicture}\n"
+                                   "      \\begin{loglogaxis}[\n"
+                                   "        width=0.6\\textwidth,\n"
+                                   "        height=0.45\\textwidth,\n"
+                                   "        xlabel={$h_{sk}$},\n"
+                                   "        xtick scale label code/.code={},\n"
+                                   "        legend pos=outer north east,\n"
+                                   "        grid style=dashed,\n"
+                                   "        cycle list name=mycycle,\n"
+                                   "      ]\n";
+
+                    ieff_graph << "    \\begin{figure}[ht!]\n"
+                                  "      \\centering\n"
+                                  "      \\begin{tikzpicture}\n"
+                                  "      \\begin{semilogxaxis}[\n"
+                                  "        width=0.6\\textwidth,\n"
+                                  "        height=0.4\\textwidth,\n"
+                                  "        xlabel={$h_{sk}$},\n"
+                                  "        %ymin=1.0, ymax=1.1,\n"
+                                  "        xtick scale label code/.code={},\n"
+                                  "        legend pos=outer north east,\n"
+                                  "        grid style=dashed,\n"
+                                  "        cycle list name=mycycle,\n"
+                                  "        y tick label style={\n"
+                                  "        /pgf/number format/.cd,\n"
+                                  "        fixed,\n"
+                                  "        fixed zerofill,\n"
+                                  "        precision=2,\n"
+                                  "        /tikz/.cd\n"
+                                  "      },\n"
+                                  "      ]\n";
+                    int current_int_ref = it->n_internal_ref;
+                    bool is_first = true;
+                    while (it->n_internal_ref == current_int_ref) {
+                        if (is_first) {
+                            estimated_error_pts << "      \\addplot\n      coordinates{\n";
+                            exact_error_pts << "      \\addplot\n      coordinates{\n";
+                            ieff_graph << "      \\addplot\n      coordinates{\n";
+                            is_first = false;
+                            estimated_error_legend << "      $h_{in} = h_{sk}/"
+                                                   << std::to_string((int)pow(2, it->n_internal_ref))
+                                                   << "$ (estimated),\n";
+                            exact_error_legend << "      $h_{in} = h_{sk}/"
+                                               << std::to_string((int)pow(2, it->n_internal_ref)) << "$ (exact),\n";
+                            ieff_error_legend << "      $h_{in} = h_{sk}/"
+                                              << std::to_string((int)pow(2, it->n_internal_ref)) << "$,\n";
+                        }
+
+                        estimated_error_pts << "        (1/" << it->n_coarse_div << ", " << it->estimated_error
+                                            << ")\n";
+                        exact_error_pts << "        (1/" << it->n_coarse_div << ", " << it->exact_error << ")\n";
+                        ieff_graph << "        (1/" << it->n_coarse_div << ", " << it->effectivity_index << ")\n";
+
+                        previous_n = it->n_order;
+                        previous_k = it->k_order;
+                        current_int_ref = it->n_internal_ref;
+                        current_problem = it->problem_name;
+                        if (it != result_vec.end()) {
+                            it++;
+                            if (it->n_order == previous_n && it->k_order == previous_k &&
+                                it->n_internal_ref != current_int_ref) {
+                                is_first = true;
+                                current_int_ref = it->n_internal_ref;
+                                estimated_error_pts << "      };\n";
+                                exact_error_pts << "      };\n";
+                                ieff_graph << "      };\n";
+                            }
+                        }
+                    }
+                    estimated_error_pts << "      };\n";
+                    exact_error_pts << "      };\n";
+                    ieff_graph << "      };\n";
+
+                    error_graph << estimated_error_pts.str();
+                    error_graph << exact_error_pts.str();
+                    error_graph << "      \\legend{\n"
+                                << estimated_error_legend.str() << exact_error_legend.str()
+                                << "      }\n"
+                                   "      \\end{loglogaxis}\n"
+                                   "      \\end{tikzpicture}\n"
+                                   "      \\caption{"
+                                << "Estimated and exact errors, $k = " << previous_k << "$, $n = " << previous_n
+                                << "$}\n    \\end{figure}\n\n";
+
+                    ieff_graph << "      \\legend{\n"
+                               << ieff_error_legend.str()
+                               << "      }\n"
+                                  "      \\end{semilogxaxis}\n"
+                                  "      \\end{tikzpicture}\n"
+                                  "      \\caption{"
+                                  "Effectivity index, $k = " << previous_k << "$, $n = " << previous_n
+                               << "$}\n    \\end{figure}\n\n";
+
+                    latex_text << "    % Error graph\n" << error_graph.str();
+                    latex_text << "    % Ieff graph\n" << ieff_graph.str();
+                    latex_text << "\\clearpage\n";
+                }
+            }
+        }
+    }
+    out << latex_text.str();
+    std::cout << latex_text.str();
+}
+
+void RunHighGradientAdaptivityProblem(const int n_divisions){
+
+    ProblemConfig config;
+    config.dimension = 2;
+    config.exact = new TLaplaceExample1;
+    config.exact.operator*().fExact = TLaplaceExample1::EArcTan;
+    config.problemname = "EArcTan";
+    config.dir_name = "MHMAdaptivity";
+    config.porder = 1;
+    config.hdivmais = 2;
+    config.materialids.insert(1);
+    config.bcmaterialids.insert(-1);
+    config.makepressurecontinuous = true;
+
+    int nCoarseRef = 7;
+    int nInternalRef = n_divisions;
 
-    constexpr auto epsilon = 0.04;
-    constexpr auto P = 1.8;
-    const auto x = coord[0];
-    const auto y = coord[1];
+    config.ndivisions = nCoarseRef;
 
-    REAL term_1 = 2 + P * cos(2 * M_PI * (x - 0.5) / epsilon);
-    REAL term_2 = 2 + P * cos(2 * M_PI * (y - 0.5) / epsilon);
+    std::string command = "mkdir -p " + config.dir_name;
+    system(command.c_str());
+
+    std::vector refLevelPerSubdomain;
+    int max_refinement = 0;
+    int adaptivity_step = 0;
+    while (max_refinement <= n_divisions) {
+        config.gmesh = CreateQuadGeoMesh(nCoarseRef, nInternalRef);
+
+        auto *mhm = new TPZMHMixedMeshControl(config.gmesh);
+        TPZManVector coarseIndexes;
+        ComputeCoarseIndices(config.gmesh, coarseIndexes);
+        bool definePartitionByCoarseIndexes = true;
+        CreateMHMCompMesh(mhm, config, definePartitionByCoarseIndexes, coarseIndexes, refLevelPerSubdomain);
+        if (max_refinement == 0) {
+            const auto n_subdomains = mhm->Coarse_to_Submesh().size();
+            refLevelPerSubdomain.resize(n_subdomains, 0);
+        }
+
+        {
+            std::string fileName = config.dir_name + "/" + config.problemname + "GMesh" + std::to_string(nCoarseRef) +
+                                   "x" + std::to_string(nCoarseRef) + std::to_string(adaptivity_step) + ".vtk";
+            std::ofstream file(fileName);
+            TPZVTKGeoMesh::PrintGMeshVTK(config.gmesh, file);
+        }
+
+        SolveMHMProblem(mhm, config);
+
+        bool postProcWithHDiv = false;
+        TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm->CMesh().operator->());
+        if (!originalMesh) DebugStop();
+        TPZMHMHDivErrorEstimator estimator(*originalMesh, mhm, postProcWithHDiv);
+        estimator.SetAdaptivityStep(adaptivity_step);
+        EstimateError(config, mhm, estimator);
+
+        auto *postprocmesh = estimator.PostProcMesh();
+        if (!postprocmesh) DebugStop();
 
-    auto perm = 1 / (term_1 * term_2);
+        MHMAdaptivity(mhm, config, postprocmesh, refLevelPerSubdomain);
+        adaptivity_step++;
 
-    for (int i = 0; i < 2; i++) {
-        res_mat(i, i) = perm;
-        res_mat(i + 2, i) = 1 / perm;
+        max_refinement = *std::max_element(refLevelPerSubdomain.begin(), refLevelPerSubdomain.end());
+        std::cout << "Max ref level: " << max_refinement << '\n';
     }
 }
 
-void PeriodicProblemForcingFunction(const TPZVec &pt, TPZVec &result) { result[0] = -1; }
+void MHMAdaptivity(TPZMHMixedMeshControl *mhm, ProblemConfig &config, TPZCompMesh *postProcMesh,
+                   std::vector &refLevelPerSubdomain) {
+
+    // Column of the flux error estimate on the element solution matrix
+    const int fluxErrorEstimateCol = 3;
+
+    TPZMultiphysicsCompMesh *cmesh = dynamic_cast(mhm->CMesh().operator->());
+    if (!cmesh) DebugStop();
+
+    TPZFMatrix &sol = postProcMesh->ElementSolution();
+    TPZSolutionMatrix &elsol = postProcMesh->ElementSolution();
+    int64_t nelem = elsol.Rows();
+
+    // Iterates through element errors to get the maximum value
+    STATE maxError = 0.;
+    for (int64_t iel = 0; iel < nelem; iel++) {
+        auto * submesh = dynamic_cast(postProcMesh->ElementVec()[iel]);
+        if (!submesh) continue;
+
+        STATE submeshError = sol(iel, fluxErrorEstimateCol);
+        if (submeshError > maxError) {
+            maxError = submeshError;
+        }
+    }
+
+    std::cout << "Max error: " << maxError << "\n";
+
+    // Refines elements which error are bigger than 30% of the maximum error
+    REAL threshold = 0.5 * maxError;
+    const auto geoToMHM = mhm->GetGeoToMHMDomain();
+    const auto interfaces = mhm->GetInterfaces();
+
+    for (int64_t iel = 0; iel < nelem; iel++) {
+        auto * submesh = dynamic_cast(postProcMesh->ElementVec()[iel]);
+        if (!submesh) continue;
+
+        STATE submeshError = sol(iel, fluxErrorEstimateCol);
+        if (submeshError > threshold) {
+            TPZGeoEl * gel = submesh->Element(0)->Reference();
+            const auto submesh_id = geoToMHM[gel->Index()];
+            refLevelPerSubdomain[submesh_id]++;
+            std::cout << "Refining submesh " << submesh_id << " which error is " << submeshError << ".\n";
+        }
+    }
+
+    for (auto i = 0; i < refLevelPerSubdomain.size(); i++) {
+        if (refLevelPerSubdomain[i] != 0) {
+            std::cout << i << ": " << refLevelPerSubdomain[i] << '\n';
+        }
+    }
+    std::cout << "Done!\n";
+}
+
+void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm, TPZMHMHDivErrorEstimator& estimator) {
+
+    std::cout << "\nError Estimation processing for MHM-Hdiv problem " << std::endl;
+
+    // Error estimation
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm->CMesh().operator->());
+    if (!originalMesh) DebugStop();
+
+    estimator.SetAnalyticSolution(config.exact);
+    estimator.PotentialReconstruction();
+
+    std::string command = "mkdir -p " + config.dir_name;
+    system(command.c_str());
+
+    TPZManVector errors;
+    TPZManVector elementerrors;
+    std::stringstream outVTK;
+    outVTK << config.dir_name << "/" << config.problemname << "-" << config.ndivisions << "-" << config.ninternalref
+           << "-Errors.vtk";
+    std::string outVTKstring = outVTK.str();
+    estimator.ComputeErrors(errors, elementerrors, outVTKstring);
+
+    {
+        std::string fileName = config.dir_name + "/" + config.problemname + "-GlobalErrors.txt";
+        std::ofstream file(fileName, std::ios::app);
+        Tools::PrintErrors(file, config, errors);
+    }
+
+    result_vec.emplace_back(config.problemname,
+                            config.ninternalref,
+                            config.ndivisions,
+                            config.porder,
+                            config.hdivmais,
+                            errors[3],
+                            errors[2],
+                            sqrt(errors[4] + errors[3]) / sqrt(errors[2])
+    );
+}
diff --git a/Projects/ErrorEstimationMHM/main_MHM.cpp b/Projects/ErrorEstimationMHM/main_MHM.cpp
index 5553e0f8..6abf2d74 100644
--- a/Projects/ErrorEstimationMHM/main_MHM.cpp
+++ b/Projects/ErrorEstimationMHM/main_MHM.cpp
@@ -2,13 +2,21 @@
 // Created by Gustavo A. Batistela on 06/07/2020.
 //
 
+#include 
 #include 
 #include 
 #include 
+#include 
+#include 
 #include 
-//#include 
+#include 
+#include 
+#include 
+#include 
 #include 
 #include 
+#include 
+#include 
 
 void RunSmoothProblem();
 void RunHighGradientProblem();
@@ -440,7 +448,7 @@ void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config) {
     TPZSSpStructMatrix strmat(cmesh.operator->());
     strmat.SetNumThreads(0 /*config.n_threads*/);
 #else
-    TPZSkylineStructMatrix strmat(cmesh.operator->());
+    TPZSkylineStructMatrix strmat(cmesh.operator->());
     strmat.SetNumThreads(0);
 #endif
 
@@ -529,23 +537,28 @@ void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfi
     int dim = control.GMesh()->Dimension();
     cmesh.SetDimModel(dim);
 
-    TPZMixedPoisson *mat = new TPZMixedPoisson(1, dim);
+    auto *mat = new TPZMixedDarcyFlow(1, dim);
 
-    TPZFMatrix K(3, 3, 0), invK(3, 3, 0);
-    K.Identity();
-    invK.Identity();
+    auto exact_lambda = [config](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+        config.exact.operator*().Exact()->Execute(loc, result, deriv);
+    };
+
+    auto ff_lambda = [config](const TPZVec &loc, TPZVec &result) {
+        config.exact.operator*().ForcingFunction()->Execute(loc, result);
+    };
 
-    mat->SetExactSol(config.exact.operator*().Exact());
-    mat->SetForcingFunction(config.exact.operator*().ForcingFunction());
-    mat->SetPermeabilityTensor(K, invK);
+    mat->SetConstantPermeability(1);
+    mat->SetExactSol(exact_lambda, 8);
+    mat->SetForcingFunction(ff_lambda, 8);
 
     cmesh.InsertMaterialObject(mat);
 
     for (auto matid : config.bcmaterialids) {
-        TPZFNMatrix<1, REAL> val1(1, 1, 0.), val2(1, 1, 0.);
+        TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+        TPZManVector val2(1, 0.);
         int bctype = 0;
-        TPZBndCond *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-        bc->TPZMaterial::SetForcingFunction(config.exact.operator*().Exact());
+        auto *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+        bc->SetForcingFunctionBC(exact_lambda);
         cmesh.InsertMaterialObject(bc);
     }
 }
diff --git a/Projects/ErrorEstimationMHM/main_RefSolValidation.cpp b/Projects/ErrorEstimationMHM/main_RefSolValidation.cpp
new file mode 100644
index 00000000..fe0bfcdc
--- /dev/null
+++ b/Projects/ErrorEstimationMHM/main_RefSolValidation.cpp
@@ -0,0 +1,393 @@
+//
+// Created by Gustavo Batistela on 25/08/21.
+//
+#include "pzgmesh.h"
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+
+void RunReferenceSolutionValidationProblem();
+
+TPZGeoMesh *CreateQuadGeoMesh(int nCoarseDiv, int nInternalRef);
+
+void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfig &config);
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes);
+
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes, const std::vector &skelsToDivide);
+
+void CreateMHMCompMeshHeteroPerm(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, int nInternalRef,
+                                 bool definePartitionByCoarseIndex, TPZManVector &mhmIndexes);
+
+void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config);
+
+void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm);
+
+int main() {
+#ifdef PZ_LOG
+    TPZLogger::InitializePZLOG();
+#endif
+    gRefDBase.InitializeAllUniformRefPatterns();
+
+    RunReferenceSolutionValidationProblem();
+
+    return 0;
+}
+
+TPZGeoMesh *CreateQuadGeoMesh(int nCoarseDiv, int nInternalRef) {
+
+    TPZManVector bcIDs(4, -1);
+    TPZGeoMesh *gmesh = Tools::CreateGeoMesh(nCoarseDiv, bcIDs);
+    gmesh->SetDimension(2);
+
+    Tools::UniformRefinement(nInternalRef, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+
+    return gmesh;
+}
+
+TPZGeoMesh *CreateCubeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &coarseIndexes) {
+
+    auto *gmesh = new TPZGeoMesh();
+    gmesh->SetDimension(3);
+    int matID = 1;
+
+    // Creates matrix with node coordinates
+    const int NodeNumber = 8;
+    REAL coordinates[NodeNumber][3] = {{0., 0., 0.}, {1., 0., 0.}, {1., 1., 0.}, {0., 1., 0.},
+                                       {0., 0., 1.}, {1., 0., 1.}, {1., 1., 1.}, {0., 1., 1.}};
+
+    // Inserts coordinates in the TPZGeoMesh object
+    for (int i = 0; i < NodeNumber; i++) {
+        int64_t nodeID = gmesh->NodeVec().AllocateNewElement();
+
+        TPZVec nodeCoord(3);
+        nodeCoord[0] = coordinates[i][0];
+        nodeCoord[1] = coordinates[i][1];
+        nodeCoord[2] = coordinates[i][2];
+
+        gmesh->NodeVec()[nodeID] = TPZGeoNode(i, nodeCoord, *gmesh);
+    }
+
+    // Creates cube element
+    TPZManVector nodeIDs(8);
+    nodeIDs[0] = 0;
+    nodeIDs[1] = 1;
+    nodeIDs[2] = 2;
+    nodeIDs[3] = 3;
+    nodeIDs[4] = 4;
+    nodeIDs[5] = 5;
+    nodeIDs[6] = 6;
+    nodeIDs[7] = 7;
+    new TPZGeoElRefPattern(nodeIDs, matID, *gmesh);
+
+    // Creates boundary faces
+    nodeIDs.Resize(4);
+    matID = -1;
+    nodeIDs[0] = 0;
+    nodeIDs[1] = 1;
+    nodeIDs[2] = 2;
+    nodeIDs[3] = 3;
+    new TPZGeoElRefPattern(nodeIDs, matID, *gmesh);
+    nodeIDs[0] = 4;
+    nodeIDs[1] = 5;
+    nodeIDs[2] = 6;
+    nodeIDs[3] = 7;
+    new TPZGeoElRefPattern(nodeIDs, matID, *gmesh);
+
+    for (int i = 0; i < 4; i++) {
+        nodeIDs[0] = (0 + i) % 4;
+        nodeIDs[1] = (1 + i) % 4;
+        nodeIDs[2] = nodeIDs[1] + 4;
+        nodeIDs[3] = nodeIDs[0] + 4;
+        new TPZGeoElRefPattern(nodeIDs, matID, *gmesh);
+    }
+
+    gmesh->BuildConnectivity();
+
+    Tools::UniformRefinement(nCoarseRef, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+
+    int64_t nElem = gmesh->NElements();
+
+    coarseIndexes.clear();
+    for (int64_t i = 0; i < nElem; i++) {
+        TPZGeoEl *gel = gmesh->Element(i);
+        if (gel->Dimension() != gmesh->Dimension() || gel->NSubElements() > 0) continue;
+        coarseIndexes.Push(i);
+    }
+
+    Tools::UniformRefinement(nInternalRef, gmesh);
+
+    return gmesh;
+}
+
+TPZGeoMesh *CreateLShapeGeoMesh(int nCoarseRef, int nInternalRef, TPZStack &mhmIndexes) {
+
+    TPZVec bcIDs(8, -1);
+    TPZGeoMesh *gmesh = Tools::CreateQuadLShapeMesh(bcIDs);
+    gmesh->SetDimension(2);
+    gmesh->BuildConnectivity();
+
+    Tools::UniformRefinement(nCoarseRef, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+
+    int64_t nElem = gmesh->NElements();
+    for (int64_t i = 0; i < nElem; i++) {
+        TPZGeoEl *gel = gmesh->Element(i);
+        if (gel->Dimension() != gmesh->Dimension() || gel->NSubElements() > 0) continue;
+        mhmIndexes.Push(i);
+    }
+
+    Tools::UniformRefinement(nInternalRef, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+
+    for (int64_t i = 0; i < mhmIndexes.size(); i++) {
+        std::cout << mhmIndexes[i] << '\n';
+    }
+    std::cout << '\n';
+
+    return gmesh;
+}
+
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes, const std::vector &skelsToDivide) {
+
+    if (definePartitionByCoarseIndex) {
+        mhm->DefinePartitionbyCoarseIndices(mhmIndexes);
+    } else {
+        mhm->DefinePartition(mhmIndexes);
+    }
+
+    // Indicate material indices to the MHM control structure
+    mhm->fMaterialIds = config.materialids;
+    mhm->fMaterialBCIds = config.bcmaterialids;
+
+    // Insert the material objects in the multiphysics mesh
+    InsertMaterialsInMHMMesh(*mhm, config);
+
+    // General approximation order settings
+    mhm->SetInternalPOrder(config.porder);
+    mhm->SetSkeletonPOrder(config.porder);
+    mhm->SetHdivmaismaisPOrder(config.hdivmais);
+
+    // Refine skeleton elements
+    for (auto skelid : skelsToDivide) {
+        mhm->DivideSkeletonElement(skelid);
+    }
+
+    mhm->DivideBoundarySkeletonElements();
+    // Creates MHM mesh
+    bool substructure = true;
+    mhm->BuildComputationalMesh(substructure);
+}
+
+void CreateMHMCompMesh(TPZMHMixedMeshControl *mhm, const ProblemConfig &config, bool definePartitionByCoarseIndex,
+                       TPZManVector &mhmIndexes) {
+
+    if (definePartitionByCoarseIndex) {
+        mhm->DefinePartitionbyCoarseIndices(mhmIndexes);
+    } else {
+        mhm->DefinePartition(mhmIndexes);
+    }
+
+    // Indicate material indices to the MHM control structure
+    mhm->fMaterialIds = config.materialids;
+    mhm->fMaterialBCIds = config.bcmaterialids;
+
+    // Insert the material objects in the multiphysics mesh
+    InsertMaterialsInMHMMesh(*mhm, config);
+
+    // General approximation order settings
+    mhm->SetInternalPOrder(config.porder);
+    mhm->SetSkeletonPOrder(config.porder);
+    mhm->SetHdivmaismaisPOrder(config.hdivmais);
+
+    // Refine skeleton elements
+    mhm->DivideSkeletonElements(0);
+    mhm->DivideBoundarySkeletonElements();
+    // Creates MHM mesh
+    bool substructure = true;
+    mhm->BuildComputationalMesh(substructure);
+}
+
+void SolveMHMProblem(TPZMHMixedMeshControl *mhm, const ProblemConfig &config) {
+
+    TPZAutoPointer cmesh = mhm->CMesh();
+
+    bool shouldrenumber = true;
+    TPZLinearAnalysis an(cmesh, shouldrenumber);
+
+#ifdef PZ_USING_MKL
+    TPZSSpStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(4);
+#else
+    TPZSkylineStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(0);
+#endif
+
+    an.SetStructuralMatrix(strmat);
+    TPZStepSolver step;
+    step.SetDirect(ELDLt);
+    an.SetSolver(step);
+    std::cout << "Assembling\n";
+    an.Assemble();
+
+    std::cout << "Solving\n";
+    an.Solve();
+    std::cout << "Finished\n";
+    an.LoadSolution(); // compute internal dofs
+
+    TPZMFSolutionTransfer transfer;
+    transfer.BuildTransferData(cmesh.operator->());
+    transfer.TransferFromMultiphysics();
+
+    TPZStack scalnames, vecnames;
+    TPZMaterial *mat = cmesh->FindMaterial(1);
+    if (!mat) {
+        DebugStop();
+    }
+
+    scalnames.Push("Pressure");
+    scalnames.Push("Permeability");
+    vecnames.Push("Flux");
+
+    if (config.exact) {
+        scalnames.Push("ExactPressure");
+        vecnames.Push("ExactFlux");
+    }
+
+    std::cout << "Post Processing...\n";
+
+    int resolution = 2;
+    std::stringstream plotname;
+    plotname << config.dir_name << "/" << config.problemname << "-" << config.ndivisions << "-" << config.ninternalref
+             << "-Results.vtk";
+    an.DefineGraphMesh(cmesh->Dimension(), scalnames, vecnames, plotname.str());
+    an.PostProcess(resolution, cmesh->Dimension());
+}
+
+void EstimateError(ProblemConfig &config, TPZMHMixedMeshControl *mhm) {
+
+    std::cout << "\nError Estimation processing for MHM-Hdiv problem " << std::endl;
+
+    // Error estimation
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm->CMesh().operator->());
+    if (!originalMesh) DebugStop();
+
+    bool postProcWithHDiv = false;
+    TPZMHMHDivErrorEstimator ErrorEstimator(*originalMesh, mhm, postProcWithHDiv);
+    ErrorEstimator.SetAnalyticSolution(config.exact);
+    ErrorEstimator.PotentialReconstruction();
+
+    std::string command = "mkdir -p " + config.dir_name;
+    system(command.c_str());
+
+    TPZManVector errors;
+    TPZManVector elementerrors;
+    std::stringstream outVTK;
+    outVTK << config.dir_name << "/" << config.problemname << "-" << config.ndivisions << "-" << config.ninternalref
+           << "-Errors.vtk";
+    std::string outVTKstring = outVTK.str();
+    ErrorEstimator.ComputeErrors(errors, elementerrors, outVTKstring);
+
+    {
+        std::string fileName = config.dir_name + "/" + config.problemname + "-GlobalErrors.txt";
+        std::ofstream file(fileName, std::ios::app);
+        Tools::PrintErrors(file, config, errors);
+    }
+
+}
+
+void InsertMaterialsInMHMMesh(TPZMHMixedMeshControl &control, const ProblemConfig &config) {
+    TPZCompMesh &cmesh = control.CMesh();
+
+    int dim = control.GMesh()->Dimension();
+    cmesh.SetDimModel(dim);
+
+    auto *mat = new TPZMixedDarcyFlow(1, dim);
+
+    auto ff_lambda = [config](const TPZVec &loc, TPZVec &result) {
+        config.exact.operator*().ForcingFunction()->Execute(loc, result);
+    };
+    auto exact_sol_lambda = [config](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+        config.exact.operator*().Exact()->Execute(loc, result, deriv);
+    };
+
+    mat->SetForcingFunction(ff_lambda, 5);
+    mat->SetExactSol(exact_sol_lambda, 5);
+    mat->SetConstantPermeability(1);
+
+    cmesh.InsertMaterialObject(mat);
+
+    for (auto matid : config.bcmaterialids) {
+        TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+        TPZManVector val2(1, 0.);
+        int bctype = 0;
+        auto *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+
+        auto ff_bc_lambda = [config](const TPZVec &loc, TPZVec &rhsVal, TPZFMatrix &matVal) {
+            config.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_bc_lambda);
+        cmesh.InsertMaterialObject(bc);
+    }
+}
+
+void RunReferenceSolutionValidationProblem() {
+
+    ProblemConfig config;
+    config.dimension = 2;
+    config.exact = new TLaplaceExample1;
+    config.exact.operator*().fExact = TLaplaceExample1::EArcTan;
+    config.problemname = "EArcTan";
+    config.dir_name = "MHMAdaptivity";
+    config.porder = 1;
+    config.hdivmais = 2;
+    config.materialids.insert(1);
+    config.bcmaterialids.insert(-1);
+    config.makepressurecontinuous = true;
+
+    int nCoarseRef = 4;
+    int nInternalRef = 4;
+
+    config.ndivisions = nCoarseRef;
+
+    std::string command = "mkdir -p " + config.dir_name;
+    system(command.c_str());
+
+    std::vector skelsToRefine;
+
+    config.gmesh = CreateQuadGeoMesh(nCoarseRef, nInternalRef);
+
+    auto *mhm = new TPZMHMixedMeshControl(config.gmesh);
+    TPZManVector coarseIndexes;
+    ComputeCoarseIndices(config.gmesh, coarseIndexes);
+    bool definePartitionByCoarseIndexes = true;
+    CreateMHMCompMesh(mhm, config, definePartitionByCoarseIndexes, coarseIndexes, skelsToRefine);
+
+    {
+        std::string fileName = config.dir_name + "/" + config.problemname + "GMesh.vtk";
+        std::ofstream file(fileName);
+        TPZVTKGeoMesh::PrintGMeshVTK(config.gmesh, file);
+    }
+
+    SolveMHMProblem(mhm, config);
+
+    bool postProcWithHDiv = false;
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm->CMesh().operator->());
+    if (!originalMesh) DebugStop();
+    TPZMHMHDivErrorEstimator estimator(*originalMesh, mhm, postProcWithHDiv);
+    EstimateError(config, mhm);
+}
diff --git a/Projects/SPE10/CMakeLists.txt b/Projects/SPE10/CMakeLists.txt
index 77ab7635..d20d4c76 100644
--- a/Projects/SPE10/CMakeLists.txt
+++ b/Projects/SPE10/CMakeLists.txt
@@ -1,7 +1,7 @@
-add_executable(SPE10 mainSPE10.cpp)
-target_link_libraries(SPE10 PRIVATE Tools ErrorEstimationLib)
+add_executable(SPE10Adaptivity mainSPE10.cpp ToolsSPE10.cpp ToolsSPE10.h TPZMultiscaleGridGen2D.cpp TPZMultiscaleGridGen2D.h)
+target_link_libraries(SPE10Adaptivity PRIVATE Tools ErrorEstimationLib)
 
-find_package(libInterpolate REQUIRED)
-target_link_libraries(SPE10 PRIVATE libInterpolate::Interpolate)
+add_executable(SPE10RefSol mainSPE10RefSol.cpp ToolsSPE10.cpp ToolsSPE10.h TPZMultiscaleGridGen2D.cpp TPZMultiscaleGridGen2D.h)
+target_link_libraries(SPE10RefSol PRIVATE Tools ErrorEstimationLib)
 
 configure_file(InputData/spe_perm.dat InputData/spe_perm.dat COPYONLY)
diff --git a/Projects/SPE10/TPZMultiscaleGridGen2D.cpp b/Projects/SPE10/TPZMultiscaleGridGen2D.cpp
new file mode 100644
index 00000000..5a96820a
--- /dev/null
+++ b/Projects/SPE10/TPZMultiscaleGridGen2D.cpp
@@ -0,0 +1,225 @@
+//
+// Created by Gustavo Batistela on 11/8/21.
+//
+
+#include "TPZMultiscaleGridGen2D.h"
+#include 
+#include 
+
+TPZMultiscaleGridGen2D::TPZMultiscaleGridGen2D(const TPZVec &minX, const TPZVec &maxX,
+                                               const TPZVec &NDivFineGrid, int NElemCoarseGrid)
+    : fMinX(minX), fMaxX(maxX), fNDivFineGrid(NDivFineGrid), fNElemCoarseGrid(NElemCoarseGrid) {
+
+    if (NDivFineGrid[0] < NElemCoarseGrid || NDivFineGrid[1] < NElemCoarseGrid) DebugStop();
+
+    fRefTreeDesiredSize = new RefTree(NElemCoarseGrid);
+
+    std::div_t div_x = std::div(NDivFineGrid[0], NElemCoarseGrid);
+    std::div_t div_y = std::div(NDivFineGrid[1], NElemCoarseGrid);
+    if (div_x.rem != 0) fRefTreeRemainderX = new RefTree(div_x.rem);
+    if (div_y.rem != 0) fRefTreeRemainderY = new RefTree(div_y.rem);
+
+    GenerateRefPatterns();
+    CreateFineGridMesh();
+    CreateSkeletonElements();
+    RefineSkeletonElements();
+    SwapSkeletonNodes();
+}
+
+TPZRefPattern TPZMultiscaleGridGen2D::CreateNonUniformLineRefPattern(const int a, const int b) {
+
+    auto *gmesh = new TPZGeoMesh();
+    gmesh->SetDimension(1);
+
+    TPZManVector coord(3, 0.);
+
+    gmesh->NodeVec().Resize(3);
+    gmesh->NodeVec()[0].Initialize(coord, *gmesh);
+
+    coord[0] = a;
+    gmesh->NodeVec()[1].Initialize(coord, *gmesh);
+    coord[0] = a + b;
+    gmesh->NodeVec()[2].Initialize(coord, *gmesh);
+
+    constexpr int mat_id = 1;
+
+    // Inserts line elements
+    TPZManVector nodes_id_vec(2, 0);
+    nodes_id_vec[1] = 2;
+    new TPZGeoElRefPattern(nodes_id_vec, mat_id, *gmesh);
+
+    nodes_id_vec[0] = 0;
+    nodes_id_vec[1] = 1;
+    new TPZGeoElRefPattern(nodes_id_vec, mat_id, *gmesh);
+    nodes_id_vec[0] = 1;
+    nodes_id_vec[1] = 2;
+    new TPZGeoElRefPattern(nodes_id_vec, mat_id, *gmesh);
+
+    gmesh->BuildConnectivity();
+    return TPZRefPattern(*gmesh);
+}
+
+void TPZMultiscaleGridGen2D::GenerateRefPatterns() {
+    std::set> ref_levels;
+
+    std::function VisitNode = [&](RefTree *node) -> void {
+        if (node->fActualSize != 1) {
+            const auto left_size = node->fChildLeft->fActualSize;
+            const auto right_size = node->fChildRight->fActualSize;
+            if (left_size != right_size) {
+                ref_levels.insert({node->fChildLeft->fActualSize, node->fChildRight->fActualSize});
+            } else {
+                ref_levels.insert({1, 1});
+            }
+            VisitNode(node->fChildLeft);
+            VisitNode(node->fChildRight);
+        }
+    };
+
+    VisitNode(fRefTreeDesiredSize);
+    if (fRefTreeRemainderX) VisitNode(fRefTreeRemainderX);
+    if (fRefTreeRemainderY) VisitNode(fRefTreeRemainderY);
+
+    for (auto it : ref_levels) {
+        const auto a = it.first;
+        const auto b = it.second;
+
+        TPZRefPattern refPattern = CreateNonUniformLineRefPattern(a, b);
+        fRefPatterns.insert({{it.first, it.second}, refPattern});
+    }
+}
+
+void TPZMultiscaleGridGen2D::CreateFineGridMesh() {
+
+    const TPZManVector bcIDs = {-2, -1, -2, -2};
+
+    TPZGenGrid2D gen(fNDivFineGrid, fMinX, fMaxX, 1, 0);
+    gen.SetRefpatternElements(true);
+
+    fGeoMesh = new TPZGeoMesh;
+    gen.Read(fGeoMesh);
+
+    gen.SetBC(fGeoMesh, 4, bcIDs[0]);
+    gen.SetBC(fGeoMesh, 5, bcIDs[1]);
+    gen.SetBC(fGeoMesh, 6, bcIDs[2]);
+    gen.SetBC(fGeoMesh, 7, bcIDs[3]);
+
+    fGeoMesh->SetDimension(2);
+}
+
+void TPZMultiscaleGridGen2D::CreateSkeletonElements() {
+
+    const auto node_vec = fGeoMesh->NodeVec();
+
+    const std::div_t div_x = std::div(fNDivFineGrid[0], fNElemCoarseGrid);
+    const std::div_t div_y = std::div(fNDivFineGrid[1], fNElemCoarseGrid);
+
+    const int n_full_size_x = div_x.quot;
+    const int n_full_size_y = div_y.quot;
+
+    const int ny_correction = div_y.rem == 0 ? 1 : 0;
+    const int nx_correction = div_x.rem == 0 ? 1 : 0;
+
+    const int n_nodes_x = fNDivFineGrid[0] + 1;
+
+    TPZManVector node_id_vec(2, 0);
+    int node0_id, node1_id;
+    for (auto iy = 0; iy < n_full_size_y - ny_correction; iy++) {
+        for (auto ix = 0; ix < n_full_size_x; ix++) {
+            node0_id =
+                fNElemCoarseGrid * fNDivFineGrid[0] + (ix + 1) * fNElemCoarseGrid + fNElemCoarseGrid * n_nodes_x * iy;
+            node1_id = node0_id + fNElemCoarseGrid;
+            node_id_vec[0] = node0_id;
+            node_id_vec[1] = node1_id;
+            auto * gel = new TPZGeoElRefPattern(node_id_vec, 666, *fGeoMesh);
+            fSkelIdToRefTree.insert({gel->Id(), fRefTreeDesiredSize});
+        }
+        if (div_x.rem != 0) {
+            node0_id = fNElemCoarseGrid * fNDivFineGrid[0] + (n_full_size_x + 1) * fNElemCoarseGrid +
+                    fNElemCoarseGrid * n_nodes_x * iy;
+            node1_id = node0_id + div_x.rem;
+            node_id_vec[0] = node0_id;
+            node_id_vec[1] = node1_id;
+            auto * gel = new TPZGeoElRefPattern(node_id_vec, 666, *fGeoMesh);
+            fSkelIdToRefTree.insert({gel->Id(), fRefTreeRemainderX});
+        }
+    }
+
+    for (auto ix = 0; ix < n_full_size_x - nx_correction; ix++) {
+        for (auto iy = 0; iy < n_full_size_y; iy++) {
+            node0_id = (fNElemCoarseGrid) * (ix + 1) + n_nodes_x * iy * fNElemCoarseGrid;
+            node1_id = node0_id + n_nodes_x * fNElemCoarseGrid;
+            node_id_vec[0] = node0_id;
+            node_id_vec[1] = node1_id;
+            auto * gel = new TPZGeoElRefPattern(node_id_vec, 666, *fGeoMesh);
+            fSkelIdToRefTree.insert({gel->Id(), fRefTreeDesiredSize});
+        }
+        if (div_y.rem != 0) {
+            node0_id = (fNElemCoarseGrid) * (ix + 1) + n_nodes_x * n_full_size_y * fNElemCoarseGrid;
+            node1_id = node0_id + n_nodes_x * div_y.rem;
+            node_id_vec[0] = node0_id;
+            node_id_vec[1] = node1_id;
+            auto * gel = new TPZGeoElRefPattern(node_id_vec, 666, *fGeoMesh);
+            fSkelIdToRefTree.insert({gel->Id(), fRefTreeRemainderY});
+        }
+    }
+}
+
+void TPZMultiscaleGridGen2D::RefineSkeletonElements() {
+
+    std::function)> RefineSkeleton = [&](std::pair data) -> void {
+        auto * node = data.second;
+        if (node->fActualSize != 1) {
+
+            auto right_size = node->fChildRight->fActualSize;
+            auto left_size = node->fChildLeft->fActualSize;
+
+            if (right_size == left_size) {
+                right_size = 1;
+                left_size = 1;
+            }
+            const auto ref_pattern = fRefPatterns.find({left_size, right_size});
+            if (ref_pattern == fRefPatterns.end()) DebugStop();
+            TPZAutoPointer ref_pattern_ptr(&ref_pattern->second);
+            auto * skel = fGeoMesh->Element(data.first);
+            skel->SetRefPattern(ref_pattern_ptr);
+            TPZManVector sons(2, nullptr);
+            skel->Divide(sons);
+
+            RefineSkeleton({sons[0]->Id(), node->fChildLeft});
+            RefineSkeleton({sons[1]->Id(), node->fChildRight});
+        }
+    };
+
+    std::for_each(fSkelIdToRefTree.begin(), fSkelIdToRefTree.end(), RefineSkeleton);
+}
+
+void TPZMultiscaleGridGen2D::SwapSkeletonNodes() {
+
+    const auto x_step = (fMaxX[0] - fMinX[0]) / fNDivFineGrid[0];
+    const auto y_step = (fMaxX[1] - fMinX[1]) / fNDivFineGrid[1];
+
+    auto CoordToNodeId = [&](const auto &coord) {
+        return std::lround(coord[0] / x_step) + (fNDivFineGrid[0] + 1) * std::lround(coord[1] / y_step);
+    };
+
+    for (auto i = 0; i < fGeoMesh->NElements(); i++) {
+        const auto * gel = fGeoMesh->Element(i);
+        if (!gel) continue;
+        if (gel->MaterialId() != fSkeletonMatId) continue;
+        if (!gel->Father()) continue;
+        TPZManVector coord0(3), coord1(3);
+        auto node0 = gel->Node(0);
+        auto node1 = gel->Node(1);
+
+        node0.GetCoordinates(coord0);
+        node1.GetCoordinates(coord1);
+
+        auto node0_id = CoordToNodeId(coord0);
+        auto node1_id = CoordToNodeId(coord1);
+
+        //std::cout << "Old nodes: " << node0.Id() << ", " << node1.Id() << '\n';
+        //std::cout << "New nodes: " << node0_id << ", " << node1_id << '\n';
+
+    }
+}
diff --git a/Projects/SPE10/TPZMultiscaleGridGen2D.h b/Projects/SPE10/TPZMultiscaleGridGen2D.h
new file mode 100644
index 00000000..76c1ed6d
--- /dev/null
+++ b/Projects/SPE10/TPZMultiscaleGridGen2D.h
@@ -0,0 +1,73 @@
+//
+// Created by Gustavo Batistela on 11/8/21.
+//
+
+#ifndef TPZMULTISCALEGRIDGEN2D_H
+#define TPZMULTISCALEGRIDGEN2D_H
+
+#include 
+
+struct RefTree {
+    int fActualSize{0};
+    RefTree *fChildLeft = nullptr;
+    RefTree *fChildRight = nullptr;
+
+    RefTree() = default;
+
+    explicit RefTree(const int actual_number) {
+        fActualSize = actual_number;
+        FillRefTree();
+    }
+
+    void FillRefTree() {
+        std::div_t division = std::div(fActualSize, 2);
+        if (fActualSize != 1) {
+            fChildLeft = new RefTree(division.quot);
+            fChildRight = new RefTree(division.quot + division.rem);
+        }
+    }
+};
+
+class TPZMultiscaleGridGen2D {
+
+public:
+    // Constructors
+    TPZMultiscaleGridGen2D() = delete;
+
+    TPZMultiscaleGridGen2D(const TPZVec &minX, const TPZVec &maxX, const TPZVec &NDivFineGrid,
+                           int NElemCoarseGrid);
+
+    [[nodiscard]] auto* GeoMesh() const { return fGeoMesh; }
+
+private:
+    // Member variables
+    const int fSkeletonMatId{-99};
+
+    const TPZManVector fNDivFineGrid{0};
+    const int fNElemCoarseGrid{0};
+    const TPZManVector fMinX{0};
+    const TPZManVector fMaxX{0};
+
+    RefTree *fRefTreeDesiredSize = nullptr;
+    RefTree *fRefTreeRemainderX = nullptr;
+    RefTree *fRefTreeRemainderY = nullptr;
+    std::map, TPZRefPattern> fRefPatterns;
+    std::map fSkelIdToRefTree;
+
+    TPZGeoMesh *fGeoMesh = nullptr;
+
+    // Private member functions
+    static TPZRefPattern CreateNonUniformLineRefPattern(int a, int b);
+
+    void GenerateRefPatterns();
+
+    void CreateFineGridMesh();
+
+    void CreateSkeletonElements();
+
+    void RefineSkeletonElements();
+
+    void SwapSkeletonNodes();
+};
+
+#endif // TPZMULTISCALEGRIDGEN2D_H
diff --git a/Projects/SPE10/ToolsSPE10.cpp b/Projects/SPE10/ToolsSPE10.cpp
new file mode 100644
index 00000000..48b3a0c6
--- /dev/null
+++ b/Projects/SPE10/ToolsSPE10.cpp
@@ -0,0 +1,422 @@
+//
+// Created by Gustavo Batistela on 10/26/21.
+//
+
+#include "Tools.h"
+#include "ToolsMHM.h"
+#include "ToolsSPE10.h"
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+
+[[maybe_unused]] TPZGeoMesh *SPE10::CreateFineGridGeoMesh() {
+
+    TPZManVector bcIDs = {-2, -1, -2, -2};
+    TPZManVector n_elements = {220, 60};
+    TPZManVector x0 = {0, 0, 0};
+    TPZManVector x1 = {220, 60, 0};
+
+    TPZGenGrid2D gen(n_elements, x0, x1, 1, 0);
+    gen.SetRefpatternElements(true);
+
+    auto *gmesh = new TPZGeoMesh;
+    gen.Read(gmesh);
+
+    gen.SetBC(gmesh, 4, bcIDs[0]);
+    gen.SetBC(gmesh, 5, bcIDs[1]);
+    gen.SetBC(gmesh, 6, bcIDs[2]);
+    gen.SetBC(gmesh, 7, bcIDs[3]);
+
+    gmesh->SetDimension(2);
+
+    return gmesh;
+}
+
+[[maybe_unused]] TPZGeoMesh *SPE10::CreateMHMGeoMesh() {
+
+    auto *gmesh = new TPZGeoMesh();
+    gmesh->SetDimension(2);
+
+    auto *gmesh2x1 = CreateRefinementGeoMesh(2, 1);
+    TPZRefPattern ref_pat2x1(*gmesh2x1);
+    TPZAutoPointer ref2x1(&ref_pat2x1);
+
+    auto *gmesh1x2 = CreateRefinementGeoMesh(1, 2);
+    TPZRefPattern ref_pat1x2(*gmesh1x2);
+    TPZAutoPointer ref1x2(&ref_pat1x2);
+
+    auto *gmesh1x1 = CreateRefinementGeoMesh(1, 1);
+    TPZRefPattern ref_pat1x1(*gmesh1x1);
+    TPZAutoPointer ref1x1(&ref_pat1x1);
+
+    TPZManVector coord(3, 0.);
+    for (int y = 0; y <= 8; y++) {
+        for (int x = 0; x <= 28; x++) {
+            coord = {8.0 * x, 8.0 * y, 0};
+            if (x == 28) coord[0] = 220.;
+            if (y == 8) coord[1] = 60.;
+
+            // Create new node
+            const auto newID = gmesh->NodeVec().AllocateNewElement();
+            gmesh->NodeVec()[newID].Initialize(coord, *gmesh);
+        }
+    }
+
+    constexpr int porousMediaMatId = 1;
+    constexpr int bcMatId1 = -1;
+    constexpr int bcMatId2 = -2;
+
+    // Inserts quad elements
+    TPZManVector nodesIdVec(4);
+    for (int y = 0; y < 8; y++) {
+        for (int x = 0; x < 28; x++) {
+            nodesIdVec[0] = 29 * y + x;
+            nodesIdVec[1] = 29 * y + x + 1;
+            nodesIdVec[2] = 29 * (y + 1) + x + 1;
+            nodesIdVec[3] = 29 * (y + 1) + x + 0;
+            auto *gel = new TPZGeoElRefPattern(nodesIdVec, porousMediaMatId, *gmesh);
+            if (x == 27 && y != 7) gel->SetRefPattern(ref1x2);
+            if (x != 27 && y == 7) gel->SetRefPattern(ref2x1);
+            if (x == 27 && y == 7) gel->SetRefPattern(ref1x1);
+        }
+    }
+
+    auto *gmesh1x1line = CreateLineRefinementGeoMesh(1);
+    TPZRefPattern ref_pat1x1line(*gmesh1x1line);
+    TPZAutoPointer ref1x1line(&ref_pat1x1line);
+
+    nodesIdVec.resize(2);
+    for (int x = 0; x < 28; x++) {
+        nodesIdVec[0] = x;
+        nodesIdVec[1] = x + 1;
+        auto *bc1 = new TPZGeoElRefPattern(nodesIdVec, bcMatId1, *gmesh);
+        nodesIdVec[0] = x + 8 * 29;
+        nodesIdVec[1] = x + 8 * 29 + 1;
+        auto *bc2 = new TPZGeoElRefPattern(nodesIdVec, bcMatId1, *gmesh);
+        if (x == 27) {
+            bc1->SetRefPattern(ref1x1line);
+            bc2->SetRefPattern(ref1x1line);
+        }
+    }
+
+    for (int y = 0; y < 8; y++) {
+        nodesIdVec[0] = y * 29;
+        nodesIdVec[1] = (y + 1) * 29;
+        auto *bc1 = new TPZGeoElRefPattern(nodesIdVec, bcMatId1, *gmesh);
+        nodesIdVec[0] = y * 29 + 28;
+        nodesIdVec[1] = (y + 1) * 29 + 28;
+        auto *bc2 = new TPZGeoElRefPattern(nodesIdVec, bcMatId2, *gmesh);
+        if (y == 7) {
+            bc1->SetRefPattern(ref1x1line);
+            bc2->SetRefPattern(ref1x1line);
+        }
+    }
+
+    gmesh->BuildConnectivity();
+
+    TPZManVector sons;
+    for (int div = 0; div < 3; div++) {
+        auto nelem = gmesh->NElements();
+        for (int64_t i = 0; i < nelem; i++) {
+            auto *gel = gmesh->ElementVec()[i];
+            const int has_sub = gel->HasSubElement();
+            if (has_sub == 0) {
+                gel->Divide(sons);
+            }
+        }
+    }
+    return gmesh;
+}
+
+[[maybe_unused]] TPZGeoMesh *SPE10::CreateRefinementGeoMesh(const int nx, const int ny) {
+    auto *gmesh = new TPZGeoMesh();
+    gmesh->SetDimension(2);
+    TPZManVector coord(3, 0.);
+
+    for (int y = 0; y <= ny; y++) {
+        for (int x = 0; x <= nx; x++) {
+            coord = {static_cast(x), static_cast(y), 0.};
+            // Create new node
+            const auto newID = gmesh->NodeVec().AllocateNewElement();
+            gmesh->NodeVec()[newID].Initialize(coord, *gmesh);
+        }
+    }
+
+    constexpr int matId = 1;
+
+    // Inserts quad elements
+    TPZManVector nodesIdVec(4, 0);
+    nodesIdVec[1] = nx;
+    nodesIdVec[2] = (ny + 1) * (nx + 1) - 1;
+    nodesIdVec[3] = ny * (nx + 1);
+    auto *father_gel = new TPZGeoElRefPattern(nodesIdVec, matId, *gmesh);
+
+    for (int y = 0; y < ny; y++) {
+        for (int x = 0; x < nx; x++) {
+            nodesIdVec[0] = (nx + 1) * y + x;
+            nodesIdVec[1] = (nx + 1) * y + x + 1;
+            nodesIdVec[2] = (nx + 1) * (y + 1) + x + 1;
+            nodesIdVec[3] = (nx + 1) * (y + 1) + x + 0;
+            auto *gel = new TPZGeoElRefPattern(nodesIdVec, matId, *gmesh);
+            gel->SetFather(father_gel);
+        }
+    }
+
+    gmesh->BuildConnectivity();
+    return gmesh;
+}
+
+[[maybe_unused]] TPZGeoMesh *SPE10::CreateLineRefinementGeoMesh(const int nx) {
+    auto *gmesh = new TPZGeoMesh();
+    gmesh->SetDimension(1);
+
+    TPZManVector coord(3, 0.);
+    for (int x = 0; x <= nx; x++) {
+        coord = {static_cast(x), 0, 0.};
+        // Create new node
+        const auto newID = gmesh->NodeVec().AllocateNewElement();
+        gmesh->NodeVec()[newID].Initialize(coord, *gmesh);
+    }
+
+    constexpr int matId = 1;
+
+    // Inserts line elements
+    TPZManVector nodesIdVec(2, 0);
+    nodesIdVec[1] = nx;
+    auto *father_gel = new TPZGeoElRefPattern(nodesIdVec, matId, *gmesh);
+
+    for (int x = 0; x < nx; x++) {
+        nodesIdVec[0] = x;
+        nodesIdVec[1] = x + 1;
+        auto *gel = new TPZGeoElRefPattern(nodesIdVec, matId, *gmesh);
+        gel->SetFather(father_gel);
+    }
+
+    gmesh->BuildConnectivity();
+    return gmesh;
+}
+
+STATE SPE10::PermeabilityFunction(const TPZVec &x) {
+    auto rounded_x = static_cast(x[0]);
+    auto rounded_y = static_cast(x[1]);
+    if (rounded_x == 220) rounded_x = 219;
+    if (rounded_y == 60) rounded_y = 59;
+    return perm_vec->operator[](rounded_x * 60 + rounded_y);
+}
+
+void SPE10::ReadSPE10CellPermeabilities() {
+
+    std::cout << "Reading permeability data...\n";
+
+    std::ifstream perm_file("InputData/spe_perm.dat", std::ios::in);
+    if (!perm_file) {
+        std::cerr << "Unable to open input file\n";
+        DebugStop();
+    }
+
+    perm_vec = new TPZManVector(n_cells, 0);
+
+    int cell_id = 0;
+    const auto start_line = 1 + n_cells * (layer - 1) / 6;
+
+    int line_num = 0;
+    int line_num2 = 0;
+    while (perm_file) {
+        line_num++;
+        line_num2++;
+        std::string line;
+        std::getline(perm_file, line, '\n');
+
+        if (line_num < start_line) continue;
+
+        std::stringstream stream(line);
+        for (int i = 0; i < 6; i++) {
+            stream >> perm_vec->operator[](cell_id);
+            cell_id++;
+        }
+        if (cell_id == n_cells) break;
+    }
+    std::cout << "Finished reading permeability data from input file!\n";
+}
+
+void SPE10::InsertMaterials(TPZCompMesh *cmesh) {
+
+    auto *mix = new TPZMixedDarcyFlow(1, cmesh->Dimension());
+    std::function &coord)> func = PermeabilityFunction;
+    mix->SetPermeabilityFunction(func);
+
+    TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+    TPZManVector val2(1, 0.);
+    constexpr int dirichlet_bc = 0;
+    constexpr int neumann_bc = 1;
+
+    // Pressure at reservoir boundary
+    val2[0] = 10;
+    TPZBndCond *pressure_left = mix->CreateBC(mix, -1, dirichlet_bc, val1, val2);
+    val2[0] = 0;
+    TPZBndCond *pressure_general = mix->CreateBC(mix, -2, neumann_bc, val1, val2);
+
+    cmesh->InsertMaterialObject(mix);
+    cmesh->InsertMaterialObject(pressure_left);
+    cmesh->InsertMaterialObject(pressure_general);
+}
+
+void SPE10::EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator) {
+
+    std::cout << "\nError Estimation processing for MHM-Hdiv problem " << std::endl;
+
+    // Error estimation
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm.CMesh().operator->());
+    if (!originalMesh) DebugStop();
+
+    estimator.PotentialReconstruction();
+
+    std::string command = "mkdir SPE10";
+    system(command.c_str());
+
+    TPZManVector errors;
+    TPZManVector elementerrors;
+    std::stringstream plotname;
+    plotname << "CroppedSPE10-Errors-Step" << estimator.AdaptivityStep() << ".vtk";
+    auto plotname_str = plotname.str();
+    estimator.ComputeErrors(errors, elementerrors, plotname_str);
+    std::cout << "Finished computing errors!\n";
+}
+
+void SPE10::EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator, TPZMultiphysicsCompMesh * ref_sol) {
+
+    std::cout << "\nError Estimation processing for MHM-Hdiv problem " << std::endl;
+
+    // Error estimation
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm.CMesh().operator->());
+    if (!originalMesh) DebugStop();
+
+    estimator.SetReferenceSolution(true);
+    estimator.SetReferenceSolutionMeshes(ref_sol->MeshVector()[1], ref_sol->MeshVector()[1]);
+
+    estimator.PotentialReconstruction();
+
+
+    std::string command = "mkdir SPE10";
+    system(command.c_str());
+
+    TPZManVector errors;
+    TPZManVector elementerrors;
+    std::stringstream plotname;
+    plotname << "CroppedSPE10-Errors-Step" << estimator.AdaptivityStep() << ".vtk";
+    auto plotname_str = plotname.str();
+    estimator.ComputeErrors(errors, elementerrors, plotname_str);
+    std::cout << "Finished computing errors!\n";
+}
+
+void SPE10::SolveMHMProblem(TPZMHMixedMeshControl &mhm, const int adaptivity_step) {
+
+    TPZAutoPointer cmesh = mhm.CMesh();
+
+    bool should_renumber = true;
+    TPZLinearAnalysis an(cmesh, should_renumber);
+
+#ifdef PZ_USING_MKL
+    TPZSSpStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(0);
+#else
+    TPZSkylineStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(8);
+#endif
+
+    an.SetStructuralMatrix(strmat);
+    TPZStepSolver step;
+    step.SetDirect(ELDLt);
+    an.SetSolver(step);
+
+    std::cout << "Assembling\n";
+    an.Assemble();
+
+    std::cout << "Solving\n";
+    an.Solve();
+    std::cout << "Finished\n";
+    an.LoadSolution();
+
+    TPZMFSolutionTransfer transfer;
+    transfer.BuildTransferData(cmesh.operator->());
+    transfer.TransferFromMultiphysics();
+
+    TPZStack scalnames, vecnames;
+    TPZMaterial *mat = cmesh->FindMaterial(1);
+    if (!mat) {
+        DebugStop();
+    }
+
+    scalnames.Push("Pressure");
+    scalnames.Push("Permeability");
+    vecnames.Push("Flux");
+
+    int resolution = 0;
+    std::stringstream plotname;
+    plotname << "CroppedSPE10-Results-Step" << adaptivity_step << ".vtk";
+    an.DefineGraphMesh(cmesh->Dimension(), scalnames, vecnames, plotname.str());
+    an.PostProcess(resolution, cmesh->Dimension());
+}
+
+TPZGeoMesh *SPE10::CreateSPE10CoarseGeoMesh() {
+    std::cout << "Creating SPE10 initial grid...\n";
+
+    const TPZManVector x0 = {0, 0, 0};
+    const TPZManVector x1 = {208., 48., 0.};
+    const TPZManVector ndiv = {13, 3, 0};
+
+    TPZGenGrid2D gen(ndiv, x0, x1);
+
+    gen.SetRefpatternElements(true);
+    auto gmesh = new TPZGeoMesh;
+    gen.Read(gmesh);
+
+    gen.SetBC(gmesh, 5, -1);
+    gen.SetBC(gmesh, 6, -2);
+    gen.SetBC(gmesh, 7, -2);
+    // gen.SetBC(gmesh, 4, -2);
+
+    std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
+
+    return gmesh;
+}
+
+void SPE10::CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm, const std::vector &skelsToDivide, const int nInternalRef) {
+
+    TPZGeoMesh *gmesh = mhm.GMesh().operator->();
+    TPZManVector coarse_indexes;
+    ComputeCoarseIndices(gmesh, coarse_indexes);
+
+    Tools::UniformRefinement(nInternalRef, 2, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+
+    mhm.DefinePartitionbyCoarseIndices(coarse_indexes);
+
+    // Indicate material indices to the MHM control structure
+    mhm.fMaterialIds = {1};
+    mhm.fMaterialBCIds = {-1, -2, -3};
+
+    // Insert the material objects in the multiphysics mesh
+    TPZCompMesh *cmesh = mhm.CMesh().operator->();
+    SPE10::InsertMaterials(cmesh);
+
+    // General approximation order settings
+    mhm.SetInternalPOrder(1);
+    mhm.SetSkeletonPOrder(1);
+    mhm.SetHdivmaismaisPOrder(2);
+
+    // Refine skeleton elements
+    for (auto skelid : skelsToDivide) {
+        mhm.DivideSkeletonElement(skelid);
+    }
+    mhm.DivideBoundarySkeletonElements();
+
+    // Creates MHM mesh
+    bool substructure = true;
+    mhm.BuildComputationalMesh(substructure);
+}
diff --git a/Projects/SPE10/ToolsSPE10.h b/Projects/SPE10/ToolsSPE10.h
new file mode 100644
index 00000000..33f911bf
--- /dev/null
+++ b/Projects/SPE10/ToolsSPE10.h
@@ -0,0 +1,42 @@
+//
+// Created by Gustavo Batistela on 10/26/21.
+//
+
+#ifndef TOOLSSPE10_H
+#define TOOLSSPE10_H
+
+#include "pzgmesh.h"
+#include 
+#include 
+
+namespace SPE10 {
+
+constexpr int nx = 220;
+constexpr int ny = 60;
+constexpr int layer = 36;
+constexpr int n_cells = nx * ny;
+static TPZManVector *perm_vec;
+
+[[maybe_unused]] TPZGeoMesh *CreateFineGridGeoMesh();
+[[maybe_unused]] TPZGeoMesh *CreateMHMGeoMesh();
+[[maybe_unused]] TPZGeoMesh *CreateRefinementGeoMesh(int nx, int ny);
+[[maybe_unused]] TPZGeoMesh *CreateLineRefinementGeoMesh(int nx);
+
+STATE PermeabilityFunction(const TPZVec &x);
+
+void ReadSPE10CellPermeabilities();
+
+void InsertMaterials(TPZCompMesh *cmesh);
+
+void EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator);
+void EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator, TPZMultiphysicsCompMesh * ref_sol);
+
+void SolveMHMProblem(TPZMHMixedMeshControl &mhm, int adaptivity_step);
+
+TPZGeoMesh *CreateSPE10CoarseGeoMesh();
+
+void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm, const std::vector &skelsToDivide, int nInternalRef);
+
+} // namespace SPE10
+
+#endif // TOOLSSPE10_H
diff --git a/Projects/SPE10/mainSPE10.cpp b/Projects/SPE10/mainSPE10.cpp
index 7a3482f0..be8cc715 100644
--- a/Projects/SPE10/mainSPE10.cpp
+++ b/Projects/SPE10/mainSPE10.cpp
@@ -3,71 +3,84 @@
 //
 
 #include "Tools.h"
+#include "ToolsSPE10.h"
+#include 
 #include 
 #include 
 #include 
+#include 
+#include 
 #include 
 #include 
-#include 
 #include 
 #include 
 
-typedef _2D::BicubicInterpolator Interpolator;
-
 // Global variables
-Interpolator interpolator;
+constexpr int nx = 220;
+constexpr int ny = 60;
+constexpr int n_cells = nx * ny;
+TPZManVector perm_vec(n_cells, 1);
+
+std::map neq_error;
 
 // Function declarations
 void ReadSPE10CellPermeabilities(TPZVec*perm_vec, int layer);
-TPZGeoMesh *CreateSPE10GeoMesh();
-void PermeabilityFunction(const TPZVec &x, TPZVec &res, TPZFMatrix &res_mat);
+TPZGeoMesh *CreateSPE10CoarseGeoMesh();
+STATE PermeabilityFunction(const TPZVec &x);
 void InsertMaterials(TPZCompMesh *cmesh);
-void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm);
-void SolveMHMProblem(TPZMHMixedMeshControl &mhm);
-void EstimateError(TPZMHMixedMeshControl *mhm);
+void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm, const std::vector &skelsToDivide);
+void SolveMHMProblem(TPZMHMixedMeshControl &mhm, int adaptivity_step);
+void EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator);
+void MHMAdaptivity(TPZMHMixedMeshControl *mhm, TPZCompMesh *postProcMesh, std::vector &skelsToRefine);
 
+constexpr int adaptivity_steps = 5;
 int main() {
 
+    gRefDBase.InitializeAllUniformRefPatterns();
+
     constexpr int layer = 36;
-    constexpr int nx = 220;
-    constexpr int ny = 60;
-    constexpr int n_cells = nx * ny;
 
-    auto perm_vec = TPZManVector(n_cells, 1);
     ReadSPE10CellPermeabilities(&perm_vec, layer);
 
-    TPZGeoMesh *gmesh = CreateSPE10GeoMesh();
-    std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
-    Tools::PrintGeometry(gmesh, "SPE10GeoMesh", false, true);
-
-    std::vector x, y, perm;
-    for (int i = 0; i < nx; i++) {
-        for (int j = 0; j < ny; j++) {
-            const int cell_id = ny * i + j;
-            const double cell_perm = perm_vec[cell_id];
-            x.push_back(0.5 + i);
-            y.push_back(0.5 + j);
-            perm.push_back(cell_perm);
-        }
-    }
+    std::vector skelsToRefine;
+    for (int step = 0; step < adaptivity_steps; step++) {
+
+        TPZGeoMesh *gmesh = CreateSPE10CoarseGeoMesh();
+        std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
 
-    interpolator.setData(x.size(), x.data(), y.data(), perm.data());
+        TPZMHMixedMeshControl mhm(gmesh);
+        CreateSPE10MHMCompMesh(mhm, skelsToRefine);
+        std::stringstream filename;
+        filename << "FineCroppedSPE10GeoMesh-Step" << step;
+        Tools::PrintGeometry(gmesh, filename.str(), false, true);
 
-    TPZMHMixedMeshControl mhm(gmesh);
-    CreateSPE10MHMCompMesh(mhm);
+        SolveMHMProblem(mhm, step);
 
-    SolveMHMProblem(mhm);
-    EstimateError(&mhm);
+        bool postProcWithHDiv = false;
+        TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm.CMesh().operator->());
+        if (!originalMesh) DebugStop();
+        TPZMHMHDivErrorEstimator estimator(*originalMesh, &mhm, postProcWithHDiv);
+        estimator.SetAdaptivityStep(step);
+        EstimateError(mhm, estimator);
+
+        auto *postprocmesh = estimator.PostProcMesh();
+        if (!postprocmesh) DebugStop();
+
+        MHMAdaptivity(&mhm, postprocmesh, skelsToRefine);
+    }
+    for (auto it: neq_error) {
+        std::cout << "(" << it.first << ", " << it.second << ")\n";
+    }
 
     return 0;
 }
 
-TPZGeoMesh *CreateSPE10GeoMesh() {
+TPZGeoMesh *CreateSPE10CoarseGeoMesh() {
     std::cout << "Creating SPE10 initial grid...\n";
 
     const TPZManVector x0 = {0, 0, 0};
-    const TPZManVector x1 = {220., 60., 0.};
-    const TPZManVector ndiv = {20, 6, 0};
+    const TPZManVector x1 = {208., 48., 0.};
+    const TPZManVector ndiv = {13, 3, 0};
 
     TPZGenGrid2D gen(ndiv, x0, x1);
 
@@ -76,6 +89,9 @@ TPZGeoMesh *CreateSPE10GeoMesh() {
     gen.Read(gmesh);
 
     gen.SetBC(gmesh, 5, -1);
+    gen.SetBC(gmesh, 6, -2);
+    gen.SetBC(gmesh, 7, -2);
+    //gen.SetBC(gmesh, 4, -2);
 
     std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
 
@@ -93,8 +109,8 @@ void ReadSPE10CellPermeabilities(TPZVec *perm_vec, const int layer) {
     }
 
     int cell_id = 0;
-    const int n_cells = perm_vec->size();
-    const int start_line = 1 + n_cells * (layer - 1) / 6;
+    const auto n_cells = perm_vec->size();
+    const auto start_line = 1 + n_cells * (layer - 1) / 6;
 
     int line_num = 0;
     int line_num2 = 0;
@@ -116,32 +132,21 @@ void ReadSPE10CellPermeabilities(TPZVec *perm_vec, const int layer) {
     std::cout << "Finished reading permeability data from input file!\n";
 }
 
-void PermeabilityFunction(const TPZVec &x, TPZVec &res, TPZFMatrix &res_mat) {
-
-    for (int i = 0; i < 2; i++) {
-        auto perm = interpolator(x[0], x[1]);
-        if (perm <= 1) {
-            perm = 1;
-        } else {
-            perm += 1;
-        }
-        res_mat(i + 2, i) = 1 / perm;
-        res_mat(i, i) = perm;
-    }
-    //std::cout << "[" << x[0] << ", " << x[1] << "]\n";
-    //std::cout << "[" << res_mat(0, 0) << " " << res_mat(0, 1) << "\n";
-    //std::cout        << res_mat(1, 0) << " " << res_mat(1, 1) << "\n";
-    //std::cout        << res_mat(2, 0) << " " << res_mat(2, 1) << "\n";
-    //std::cout        << res_mat(3, 0) << " " << res_mat(3, 1) << "]\n\n";
+STATE PermeabilityFunction(const TPZVec &x) {
+    auto rounded_x = static_cast(x[0]);
+    auto rounded_y = static_cast(x[1]);
+    if (rounded_x == 220) rounded_x = 219;
+    if (rounded_y == 60) rounded_y = 59;
+    return perm_vec[rounded_x * 60 + rounded_y];
 }
 
-void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm) {
+void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm, const std::vector &skelsToDivide) {
 
     TPZGeoMesh *gmesh = mhm.GMesh().operator->();
     TPZManVector coarse_indexes;
     ComputeCoarseIndices(gmesh, coarse_indexes);
 
-    int nInternalRef = 3;
+    int nInternalRef = adaptivity_steps - 1;
     Tools::UniformRefinement(nInternalRef, 2, gmesh);
     Tools::DivideLowerDimensionalElements(gmesh);
 
@@ -158,34 +163,31 @@ void CreateSPE10MHMCompMesh(TPZMHMixedMeshControl &mhm) {
     // General approximation order settings
     mhm.SetInternalPOrder(1);
     mhm.SetSkeletonPOrder(1);
-    mhm.SetHdivmaismaisPOrder(1);
+    mhm.SetHdivmaismaisPOrder(3);
 
     // Refine skeleton elements
-    mhm.DivideSkeletonElements(0);
+    for (auto skelid : skelsToDivide) {
+        mhm.DivideSkeletonElement(skelid);
+    }
     mhm.DivideBoundarySkeletonElements();
 
     // Creates MHM mesh
     bool substructure = true;
     mhm.BuildComputationalMesh(substructure);
-    //{
-    //    std::string fileName = "CompMesh.txt";
-    //    std::ofstream file(fileName);
-    //    mhm.CMesh()->Print(file);
-    //}
 }
 
-void SolveMHMProblem(TPZMHMixedMeshControl &mhm) {
+void SolveMHMProblem(TPZMHMixedMeshControl &mhm, const int adaptivity_step) {
 
     TPZAutoPointer cmesh = mhm.CMesh();
-
+    std::cout << "*** NEQ: ," << mhm.CMesh()->NEquations() << '\n';
     bool should_renumber = true;
-    TPZAnalysis an(cmesh, should_renumber);
+    TPZLinearAnalysis an(cmesh, should_renumber);
 
 #ifdef PZ_USING_MKL
-    TPZSymetricSpStructMatrix strmat(cmesh.operator->());
-    strmat.SetNumThreads(8);
+    TPZSSpStructMatrix strmat(cmesh.operator->());
+    strmat.SetNumThreads(0);
 #else
-    TPZSkylineStructMatrix strmat(cmesh.operator->());
+    TPZSkylineStructMatrix strmat(cmesh.operator->());
     strmat.SetNumThreads(8);
 #endif
 
@@ -217,45 +219,110 @@ void SolveMHMProblem(TPZMHMixedMeshControl &mhm) {
     vecnames.Push("Flux");
 
     int resolution = 0;
-    std::string plotname = "SPE10-Results.vtk";
-    an.DefineGraphMesh(cmesh->Dimension(), scalnames, vecnames, plotname);
+    std::stringstream plotname;
+    plotname << "FineCroppedSPE10-Results-Step" << adaptivity_step << ".vtk";
+    an.DefineGraphMesh(cmesh->Dimension(), scalnames, vecnames, plotname.str());
     an.PostProcess(resolution, cmesh->Dimension());
 }
 
 void InsertMaterials(TPZCompMesh *cmesh) {
 
-    auto *mix = new TPZMixedPoisson(1, cmesh->Dimension());
-    TPZAutoPointer> perm_function = new TPZDummyFunction(&PermeabilityFunction, 3);
-    mix->SetPermeabilityFunction(perm_function);
+    auto *mix = new TPZMixedDarcyFlow(1, cmesh->Dimension());
+    std::function &coord)> func = PermeabilityFunction;
+    mix->SetPermeabilityFunction(func);
 
-    TPZFNMatrix<1, REAL> val1(1, 1, 0.), val2(1, 1, 0.);
+    TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+    TPZManVector val2(1, 0.);
     constexpr int dirichlet_bc = 0;
+    constexpr int neumann_bc = 1;
 
     // Pressure at reservoir boundary
-    val2(0, 0) = 1;
+    val2[0] = 10;
     TPZBndCond *pressure_left = mix->CreateBC(mix, -1, dirichlet_bc, val1, val2);
+    val2[0] = 0;
+    TPZBndCond *pressure_general = mix->CreateBC(mix, -2, neumann_bc, val1, val2);
 
     cmesh->InsertMaterialObject(mix);
     cmesh->InsertMaterialObject(pressure_left);
+    cmesh->InsertMaterialObject(pressure_general);
 }
 
-void EstimateError(TPZMHMixedMeshControl *mhm) {
+void EstimateError(TPZMHMixedMeshControl &mhm, TPZMHMHDivErrorEstimator &estimator) {
 
     std::cout << "\nError Estimation processing for MHM-Hdiv problem " << std::endl;
 
     // Error estimation
-    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm->CMesh().operator->());
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm.CMesh().operator->());
     if (!originalMesh) DebugStop();
 
-    bool postProcWithHDiv = true;
-    TPZMHMHDivErrorEstimator ErrorEstimator(*originalMesh, mhm, postProcWithHDiv);
-    ErrorEstimator.PotentialReconstruction();
+    estimator.PotentialReconstruction();
 
     std::string command = "mkdir SPE10";
     system(command.c_str());
 
     TPZManVector errors;
     TPZManVector elementerrors;
-    std::string outVTK = "SPE10-Errors.vtk";
-    ErrorEstimator.ComputeErrors(errors, elementerrors, outVTK);
+    std::stringstream plotname;
+    plotname << "FineCroppedSPE10-Errors-Step" << estimator.AdaptivityStep() << ".vtk";
+    auto plotname_str = plotname.str();
+    estimator.ComputeErrors(errors, elementerrors, plotname_str);
+    Tools::PrintErrors(std::cout, errors);
+
+    neq_error.insert({mhm.CMesh()->NEquations(), errors[3]});
+
+    std::cout << "Finished computing errors!\n";
 }
+
+void MHMAdaptivity(TPZMHMixedMeshControl *mhm, TPZCompMesh *postProcMesh, std::vector &skelsToRefine) {
+
+    // Column of the flux error estimate on the element solution matrix
+    const int fluxErrorEstimateCol = 3;
+
+    TPZMultiphysicsCompMesh *cmesh = dynamic_cast(mhm->CMesh().operator->());
+    if (!cmesh) DebugStop();
+
+    TPZFMatrix &sol = postProcMesh->ElementSolution();
+    TPZSolutionMatrix &elsol = postProcMesh->ElementSolution();
+    int64_t nelem = elsol.Rows();
+
+    // Iterates through element errors to get the maximum value
+    STATE maxError = 0.;
+    for (int64_t iel = 0; iel < nelem; iel++) {
+        auto * submesh = dynamic_cast(postProcMesh->ElementVec()[iel]);
+        if (!submesh) continue;
+
+        STATE submeshError = sol(iel, fluxErrorEstimateCol);
+        if (submeshError > maxError) {
+            maxError = submeshError;
+        }
+    }
+
+
+    // Refines elements which error are bigger than 30% of the maximum error
+    REAL threshold = 0.5 * maxError;
+    const auto geoToMHM = mhm->GetGeoToMHMDomain();
+    const auto interfaces = mhm->GetInterfaces();
+
+    std::set interfacesToRefine;
+    for (int64_t iel = 0; iel < nelem; iel++) {
+        auto * submesh = dynamic_cast(postProcMesh->ElementVec()[iel]);
+        if (!submesh) continue;
+
+        STATE submeshError = sol(iel, fluxErrorEstimateCol);
+        if (submeshError > threshold) {
+        //if (true) {
+            TPZGeoEl * gel = submesh->Element(0)->Reference();
+            const auto submesh_id = geoToMHM[gel->Index()];
+            //std::cout << "Refining submesh " << submesh_id << " which error is " << submeshError << ".\n";
+            for (auto interface : interfaces) {
+                if (interface.second.first == submesh_id || interface.second.second == submesh_id) {
+                    interfacesToRefine.insert(interface.first);
+                }
+            }
+        }
+    }
+
+    for (auto it : interfacesToRefine) {
+        skelsToRefine.push_back(it);
+    }
+}
\ No newline at end of file
diff --git a/Projects/SPE10/mainSPE10RefSol.cpp b/Projects/SPE10/mainSPE10RefSol.cpp
new file mode 100644
index 00000000..c3da0b0e
--- /dev/null
+++ b/Projects/SPE10/mainSPE10RefSol.cpp
@@ -0,0 +1,164 @@
+//
+// Created by Gustavo Batistela on 4/5/21.
+//
+
+#include "Tools.h"
+#include "ToolsSPE10.h"
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+
+constexpr int porder = 1;
+constexpr int korder = 2;
+constexpr int nInternalRef = 5;
+
+TPZMultiphysicsCompMesh RunFineScaleProblem();
+
+int main() {
+
+    gRefDBase.InitializeAllUniformRefPatterns();
+
+    SPE10::ReadSPE10CellPermeabilities();
+
+    auto fine_cmesh = RunFineScaleProblem();
+    // Create geomesh
+    TPZGeoMesh *gmesh = SPE10::CreateSPE10CoarseGeoMesh();
+    std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
+
+    // Create compmesh
+    TPZMHMixedMeshControl mhm(gmesh);
+
+    std::vector skelsToRefine;
+    SPE10::CreateSPE10MHMCompMesh(mhm, skelsToRefine, nInternalRef);
+    constexpr int step = 0;
+    std::stringstream filename;
+    filename << "CroppedSPE10GeoMesh-Step" << step;
+    Tools::PrintGeometry(gmesh, filename.str(), false, true);
+
+    SPE10::SolveMHMProblem(mhm, step);
+
+    // Estimate error
+    bool postProcWithHDiv = false;
+    TPZMultiphysicsCompMesh *originalMesh = dynamic_cast(mhm.CMesh().operator->());
+    if (!originalMesh) DebugStop();
+
+    TPZMHMHDivErrorEstimator estimator(*originalMesh, &mhm, postProcWithHDiv);
+    estimator.SetAdaptivityStep(step);
+    SPE10::EstimateError(mhm, estimator, &fine_cmesh);
+
+    return 0;
+}
+
+TPZMultiphysicsCompMesh RunFineScaleProblem() {
+    // Create geomesh
+    TPZGeoMesh *gmesh = SPE10::CreateSPE10CoarseGeoMesh();
+    std::cout << "SPE10 initial grid created. NElem: " << gmesh->NElements() << "\n";
+
+    Tools::UniformRefinement(nInternalRef, 2, gmesh);
+    Tools::DivideLowerDimensionalElements(gmesh);
+    const int dim = gmesh->Dimension();
+
+    TPZMultiphysicsCompMesh mixed_cmesh(gmesh);
+
+    SPE10::InsertMaterials(&mixed_cmesh);
+    mixed_cmesh.ApproxSpace().SetAllCreateFunctionsMultiphysicElem();
+
+    TPZManVector active(2, 1);
+    TPZManVector meshvector(2, 0);
+
+    auto* cmesh_flux = new TPZCompMesh(gmesh);
+    gmesh->ResetReference();
+    auto *null_mat = new TPZNullMaterial(1);
+    null_mat->SetDimension(dim);
+    cmesh_flux->InsertMaterialObject(null_mat);
+
+    for (auto matid : {-1, -2}) {
+        TPZFNMatrix<1, REAL> val1(1, 1, 0.);
+        TPZManVector val2(1, 1.);
+        const int bctype = 0;
+        auto * bc = null_mat->CreateBC(null_mat, matid, bctype, val1, val2);
+        cmesh_flux->InsertMaterialObject(bc);
+    }
+    cmesh_flux->SetDefaultOrder(porder);
+    cmesh_flux->ApproxSpace().SetAllCreateFunctionsHDiv(dim);
+    cmesh_flux->AutoBuild();
+
+    cmesh_flux->InitializeBlock();
+
+    auto *cmesh_pressure = new TPZCompMesh(gmesh);
+    auto *null_mat_p = new TPZNullMaterial(1);
+    null_mat_p->SetDimension(cmesh_pressure->Dimension());
+    cmesh_pressure->InsertMaterialObject(null_mat_p);
+
+    cmesh_pressure->SetDefaultOrder(porder + korder);
+    // cmesh_pressure->SetDefaultOrder(problem.porder);
+    cmesh_pressure->ApproxSpace().SetAllCreateFunctionsContinuous();
+    cmesh_pressure->ApproxSpace().CreateDisconnectedElements(true);
+    cmesh_pressure->AutoBuild();
+    int64_t n_connects = cmesh_pressure->NConnects();
+    for (int64_t i = 0; i < n_connects; ++i) {
+        cmesh_pressure->ConnectVec()[i].SetLagrangeMultiplier(1);
+    }
+
+    meshvector[0] = cmesh_flux;
+    meshvector[1] = cmesh_pressure;
+
+    TPZCompMeshTools::AdjustFluxPolynomialOrders(meshvector[0], korder);
+    TPZCompMeshTools::SetPressureOrders(meshvector[0], meshvector[1]);
+
+    mixed_cmesh.BuildMultiphysicsSpace(active, meshvector);
+    mixed_cmesh.LoadReferences();
+    bool keepmatrix = false;
+    bool keeponelagrangian = true;
+    TPZCompMeshTools::CreatedCondensedElements(&mixed_cmesh, keeponelagrangian, keepmatrix);
+    mixed_cmesh.InitializeBlock();
+
+    TPZHybridizeHDiv hybridizer;
+    TPZMultiphysicsCompMesh* hybridMesh = hybridizer.Hybridize(&mixed_cmesh);
+    hybridMesh->CleanUpUnconnectedNodes();
+    hybridMesh->AdjustBoundaryElements();
+
+    TPZLinearAnalysis an(hybridMesh);
+
+#ifdef PZ_USING_MKL
+    TPZSSpStructMatrix strmat(hybridMesh);
+    strmat.SetNumThreads(0);
+#else
+    TPZSkylineStructMatrix strmat(Hybridmesh);
+    strmat.SetNumThreads(0);
+#endif
+
+    std::set matIds = {1, -1, -2, hybridizer.fInterfaceMatid.first, hybridizer.fInterfaceMatid.second};
+    strmat.SetMaterialIds(matIds);
+
+    an.SetStructuralMatrix(strmat);
+
+    auto *direct = new TPZStepSolver;
+    direct->SetDirect(ELDLt);
+    an.SetSolver(*direct);
+    delete direct;
+    direct = nullptr;
+    an.Assemble();
+    an.Solve();
+
+    TPZStack scalnames, vecnames;
+    scalnames.Push("ExactPressure");
+    scalnames.Push("Pressure");
+    scalnames.Push("Permeability");
+    vecnames.Push("ExactFlux");
+    vecnames.Push("Flux");
+
+    std::stringstream sout;
+    sout << "SPE10FINE-Results.vtk";
+    an.DefineGraphMesh(2, scalnames, vecnames, sout.str());
+    constexpr int resolution = 0;
+    an.PostProcess(resolution, hybridMesh->Dimension());
+    return mixed_cmesh;
+}
diff --git a/Projects/Tools/Tools.cpp b/Projects/Tools/Tools.cpp
index 517d3645..c6779bbc 100644
--- a/Projects/Tools/Tools.cpp
+++ b/Projects/Tools/Tools.cpp
@@ -13,9 +13,12 @@
 #include 
 #include 
 #include 
+#include 
+#include 
+#include 
+#include 
 #include "DataStructure.h"
-#include "TPZNullMaterial.h"
-#include "DarcyFlow/TPZMixedDarcyFlow.h"
+
 
 #include "pzelementgroup.h"
 
@@ -36,9 +39,10 @@ void Tools::PrintGeometry(TPZGeoMesh *gmesh, const std::string &file_name, bool
 
 TPZCompMesh* Tools::CreatePressureMesh(const ProblemConfig& problem) {
     TPZCompMesh* cmesh = new TPZCompMesh(problem.gmesh);
-    TPZMaterial* mat = 0;
+    TPZNullMaterial* mat = nullptr;
     for (auto matid : problem.materialids) {
-        TPZNullMaterial<>* mix = new TPZNullMaterial<>(matid, cmesh->Dimension());
+        auto * mix = new TPZNullMaterial(matid);
+        mix->SetDimension(cmesh->Dimension());
         if (!mat) mat = mix;
         cmesh->InsertMaterialObject(mix);
     }
@@ -57,33 +61,34 @@ TPZCompMesh* Tools::CreatePressureMesh(const ProblemConfig& problem) {
         Prefinamento(cmesh, problem.ndivisions, problem.porder);
     }
     
-    
     return cmesh;
 }
 
 TPZCompMesh* Tools::CreateFluxHDivMesh(const ProblemConfig& problem) {
     int dim = problem.gmesh->Dimension();
-    TPZCompMesh* cmesh = new TPZCompMesh(problem.gmesh);
-    TPZNullMaterial<>* mat = NULL;
+    auto* cmesh = new TPZCompMesh(problem.gmesh);
+    TPZNullMaterial* mat = nullptr;
     problem.gmesh->ResetReference();
     for (auto matid : problem.materialids) {
-        TPZNullMaterial<>* mix = new TPZNullMaterial<>(matid);
+        auto * mix = new TPZNullMaterial(matid);
         mix->SetDimension(dim);
         if (!mat) mat = mix;
         cmesh->InsertMaterialObject(mix);
     }
-    if(!mat) DebugStop();
     for (auto matid : problem.bcmaterialids) {
         TPZFNMatrix<1, REAL> val1(1, 1, 0.);
-        TPZVec val2(1, 1.);
+        TPZManVector val2(1, 1.);
         int bctype;
         if (matid == -1 || matid == 2) {
             bctype = 0;
         } else {
             bctype = 1;
         }
-        TPZBndCondT* bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-//        bc->TPZMaterial::SetForcingFunction(problem.exact.operator*().Exact());
+        auto * bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+        auto ff_lambda = [problem](const TPZVec &loc, TPZVec &rhsVal, TPZFMatrix &matVal) {
+            problem.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_lambda);
         cmesh->InsertMaterialObject(bc);
     }
     cmesh->SetDefaultOrder(problem.porder);
@@ -93,25 +98,22 @@ TPZCompMesh* Tools::CreateFluxHDivMesh(const ProblemConfig& problem) {
     if (problem.prefine) {
         Prefinamento(cmesh, problem.ndivisions, problem.porder);
     }
-    
-    
+
     cmesh->InitializeBlock();
     return cmesh;
-    
 }
 
-TPZMultiphysicsCompMesh* Tools::CreateHDivMesh(const ProblemConfig& problem) {
+TPZMultiphysicsCompMesh* Tools::CreateMixedMesh(const ProblemConfig& problem) {
 
-    TPZMultiphysicsCompMesh* cmesh = new TPZMultiphysicsCompMesh(problem.gmesh);
+    auto* cmesh = new TPZMultiphysicsCompMesh(problem.gmesh);
+    TPZMixedDarcyFlow *mat = nullptr;
     TPZFMatrix K(3, 3, 0), invK(3, 3, 0);
     K.Identity();
     invK.Identity();
     
     STATE Km = problem.Km;
 
-    
     if (problem.TensorNonConst && problem.gmesh->Dimension() == 3) {
-        DebugStop();
         for (int i = 0; i < 3; i++) {
             for (int j = 0; j < 3; j++) {
                 if (i == j) {
@@ -126,19 +128,19 @@ TPZMultiphysicsCompMesh* Tools::CreateHDivMesh(const ProblemConfig& problem) {
         
     }
 
-//    K.Print(std::cout);
-//    invK.Print(std::cout);
+    for (auto matid : problem.materialids) {
+        auto *mix = new TPZMixedDarcyFlow(matid, cmesh->Dimension());
 
-    typedef TPZMixedDarcyFlow TPZMixedPoisson;
-    
-    TPZMixedPoisson* mat = NULL;
+        auto ff_lambda = [problem](const TPZVec &loc, TPZVec &result) {
+            problem.exact.operator*().ForcingFunction()->Execute(loc, result);
+        };
+        auto exact_sol_lambda = [problem](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+            problem.exact.operator*().Exact()->Execute(loc, result, deriv);
+        };
 
-    for (auto matid : problem.materialids) {
-        TPZMixedPoisson *mix = new TPZMixedPoisson(matid, cmesh->Dimension());
-        int porder = 5;
-        mix->SetForcingFunction(problem.exact.operator*().ForceFunc(),porder);
-        mix->SetExactSol(problem.exact.operator*().ExactSolution(),porder);
-        mix->SetConstantPermeability(1.);
+        mix->SetForcingFunction(ff_lambda, 5);
+        mix->SetExactSol(exact_sol_lambda, 5);
+        mix->SetConstantPermeability(K(0,0));
 
         if (!mat) mat = mix;
 
@@ -148,31 +150,33 @@ TPZMultiphysicsCompMesh* Tools::CreateHDivMesh(const ProblemConfig& problem) {
     
     for (auto matid : problem.bcmaterialids) {
         TPZFNMatrix<1, REAL> val1(1, 1, 0.);
-        TPZManVector val2(1, 0.);
+        TPZManVector val2(1, 0.);
         int bctype;
     
         switch (matid) {
-            case -1 :{
-            bctype = 0;
+            case -1 : {
+                bctype = 0;
                 break;
             }
-                
-                
-            case -2:{
-            bctype = 1;
-    
-            break;
-            }
-            case -3:{
-            bctype = 4;// different from mixed (bctype 2) already implemented on TPZMixedPoisson3d
-            val1(0,0) = Km ;
+            case -2: {
+                bctype = 1;
 
-                
-            break;
+                break;
+            }
+            case -3: {
+                bctype = 4;// different from mixed (bctype 2) already implemented on TPZMixedPoisson3d
+                val1(0, 0) = Km;
+                break;
+            }
+            default: {
+                bctype = -1;
             }
         }
-        TPZBndCondT* bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-        bc->SetForcingFunctionBC(problem.exact.operator*().ExactSolution());
+        auto * bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+        auto ff_lambda = [problem](const TPZVec &loc, TPZVec &rhsVal, TPZFMatrix &matVal) {
+            problem.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_lambda);
         cmesh->InsertMaterialObject(bc);
     }
     cmesh->ApproxSpace().SetAllCreateFunctionsMultiphysicElem();
@@ -417,14 +421,14 @@ void Tools::SolveHybridProblem(TPZCompMesh *Hybridmesh, std::pair Inte
     TPZLinearAnalysis an(Hybridmesh);
 
 #ifdef PZ_USING_MKL
-    TPZSSpStructMatrix<> strmat(Hybridmesh);
+    TPZSSpStructMatrix strmat(Hybridmesh);
     strmat.SetNumThreads(0);
     //        strmat.SetDecomposeType(ELDLt);
 #else
     //    TPZFrontStructMatrix > strmat(Hybridmesh);
     //    strmat.SetNumThreads(2);
     //    strmat.SetDecomposeType(ELDLt);
-    TPZSkylineStructMatrix strmat(Hybridmesh);
+    TPZSkylineStructMatrix strmat(Hybridmesh);
     strmat.SetNumThreads(0);
 #endif
 
@@ -464,7 +468,7 @@ void Tools::SolveHybridProblem(TPZCompMesh *Hybridmesh, std::pair Inte
         sout << problem.dir_name << "/" << "OriginalHybrid_Order_" << problem.porder << "Nref_" << problem.ndivisions
              << "NAdapStep_" << problem.adaptivityStep << ".vtk";
         an.DefineGraphMesh(2, scalnames, vecnames, sout.str());
-        int resolution = 0;
+        int resolution = 2;
         an.PostProcess(resolution, Hybridmesh->Dimension());
 
         if (problem.exact.operator*().Exact()) {
@@ -518,7 +522,7 @@ void Tools::SolveHybridProblem(TPZCompMesh *Hybridmesh, std::pair Inte
 }
 
 void Tools::SolveMixedProblem(TPZCompMesh* cmesh_HDiv, const ProblemConfig& config) {
-#ifdef ERRORESTIMATION_DEBUG
+#ifdef PZDEBUG
     {
         std::ofstream out("gmeshSolve.vtk");
         TPZVTKGeoMesh::PrintGMeshVTK(config.gmesh, out);
@@ -530,8 +534,8 @@ void Tools::SolveMixedProblem(TPZCompMesh* cmesh_HDiv, const ProblemConfig& conf
     TPZLinearAnalysis an(cmesh_HDiv, false);
 
 
-    TPZSSpStructMatrix<> strmat(cmesh_HDiv);
-    strmat.SetNumThreads(0);
+    TPZSSpStructMatrix strmat(cmesh_HDiv);
+    strmat.SetNumThreads(8);
     an.SetStructuralMatrix(strmat);
 
     std::set matids;
@@ -578,35 +582,35 @@ void Tools::SolveMixedProblem(TPZCompMesh* cmesh_HDiv, const ProblemConfig& conf
     int resolution = 2;
     an.PostProcess(resolution, dim);
 
-    if (config.exact.operator*().Exact()) {
-        TPZManVector errors(4, 0.);
-        an.SetThreadsForError(0);
-        an.SetExact(config.exact.operator*().ExactSolution());
-        an.PostProcessError(errors, false);
-
-        // Erro
-        std::ofstream myfile;
-        /*Error on MixedPoisson
-           [0] L2 for pressure
-           [1] L2 for flux
-           [2] L2 for div(flux)
-           [3] Grad pressure (Semi H1)
-           [4] Hdiv norm
-           */
-
-          // Erro
-          myfile.open("ErrorMixed.txt", std::ios::app);
-          myfile << "\n\n Error for Mixed formulation ";
-          myfile << "\n-------------------------------------------------- \n";
-          myfile << "Ndiv = " << config.ndivisions
-                 << " Order k = " << config.porder << " n "< errors(4, 0.);
+    //    an.SetThreadsForError(0);
+    //    an.SetExact(config.exact.operator*().ExactSolution());
+    //    an.PostProcessError(errors, false);
+
+    //    // Erro
+    //    std::ofstream myfile;
+    //    /*Error on MixedPoisson
+    //       [0] L2 for pressure
+    //       [1] L2 for flux
+    //       [2] L2 for div(flux)
+    //       [3] Grad pressure (Semi H1)
+    //       [4] Hdiv norm
+    //       */
+
+    //      // Erro
+    //      myfile.open("ErrorMixed.txt", std::ios::app);
+    //      myfile << "\n\n Error for Mixed formulation ";
+    //      myfile << "\n-------------------------------------------------- \n";
+    //      myfile << "Ndiv = " << config.ndivisions
+    //             << " Order k = " << config.porder << " n "<Dimension());
-        int porder = 5;
-        mix->SetExactSol(problem.exact.operator*().ExactSolution(),porder);
-        mix->SetForcingFunction(problem.exact.operator*().ForceFunc(),porder);
+        auto *mix = new TPZDarcyFlow(matid, cmesh->Dimension());
+
+
+        auto ff_lambda = [problem](const TPZVec &loc, TPZVec &result) {
+            problem.exact.operator*().ForcingFunction()->Execute(loc, result);
+        };
+        auto exact_sol_lambda = [problem](const TPZVec &loc, TPZVec &result, TPZFMatrix &deriv) {
+            problem.exact.operator*().Exact()->Execute(loc, result, deriv);
+        };
+
+        mix->SetForcingFunction(ff_lambda, 5);
+        mix->SetExactSol(exact_sol_lambda, 5);
 
         if (!mat) mat = mix;
         cmesh->InsertMaterialObject(mix);
@@ -674,10 +684,13 @@ TPZCompMesh* Tools::CMeshH1(ProblemConfig problem) {
 
     for (auto matid : problem.bcmaterialids) {
         TPZFNMatrix<1, REAL> val1(1, 1, 0.);
-        TPZManVector val2(1, 0.);
+        TPZManVector val2(1, 0.);
         int bctype = 0;
-        TPZBndCondT *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
-        bc->SetForcingFunctionBC(problem.exact.operator*().ExactSolution());
+        auto *bc = mat->CreateBC(mat, matid, bctype, val1, val2);
+        auto ff_lambda = [problem](const TPZVec &loc, TPZVec &rhsVal, TPZFMatrix &matVal) {
+            problem.exact.operator*().Exact()->Execute(loc, rhsVal, matVal);
+        };
+        bc->SetForcingFunctionBC(ff_lambda);
 
         cmesh->InsertMaterialObject(bc);
     }
@@ -709,7 +722,6 @@ void Tools::hAdaptivity(TPZCompMesh* postProcessMesh, TPZGeoMesh* gmeshToRefine,
         if (cel->Dimension() != postProcessMesh->Dimension()) continue;
         REAL elementError = elsol(iel, fluxErrorEstimateCol);
 
-
         if (elementError > maxError) {
             maxError = elementError;
         }
@@ -726,15 +738,11 @@ void Tools::hAdaptivity(TPZCompMesh* postProcessMesh, TPZGeoMesh* gmeshToRefine,
         if (cel->Dimension() != postProcessMesh->Dimension()) continue;
 
         REAL elementError = elsol(iel, fluxErrorEstimateCol);
-        //prefinement
         if (elementError > threshold) {
-
-            std::cout << "element error " << elementError << "el " << iel << "\n";
             TPZGeoEl* gel = cel->Reference();
-            int iel = gel->Id();
 
             TPZVec sons;
-            TPZGeoEl* gelToRefine = gmeshToRefine->ElementVec()[iel];
+            TPZGeoEl* gelToRefine = gmeshToRefine->ElementVec()[gel->Id()];
             if (gelToRefine && !gelToRefine->HasSubElement()) {
                 gelToRefine->Divide(sons);
 #ifdef LOG4CXX2
@@ -753,15 +761,6 @@ void Tools::hAdaptivity(TPZCompMesh* postProcessMesh, TPZGeoMesh* gmeshToRefine,
                 }
 #endif
             }
-        } else {
-            std::cout << "como refinar em p? " << "\n";
-//            TPZInterpolationSpace *sp = dynamic_cast(cel);
-//            if(!sp) continue;
-//            int level = sp->Reference()->Level();
-//            int ordem = config.porder + (config.adaptivityStep -1 ) + (level);
-//            std::cout<<"level "<< level<<" ordem "<PRefine(ordem);
-
         }
     }
     DivideLowerDimensionalElements(gmeshToRefine);
@@ -1022,6 +1021,7 @@ TPZGeoMesh* Tools::CreateGeoMesh(int nel, TPZVec& bcids, int dim, bool isOr
         return gmesh;
     }
     DebugStop(); // Dim should be 2 or 3
+    return nullptr;
 }
 
 void Tools::DrawGeoMesh(ProblemConfig &config, PreConfig &preConfig) {
@@ -1059,6 +1059,9 @@ void Tools::PrintErrors(std::ofstream& out, const ProblemConfig& config, const T
     if (config.adaptivityStep != -1) {
         ss << ", AdaptivityStep = " << config.adaptivityStep;
     }
+    if (config.cmesh) {
+        ss << ", NEquations = " << config.cmesh->NEquations();
+    }
     ss << '\n';
     ss << "Global estimator = " << error_vec[3] << "\n";
     ss << "|ufem-urec| = " << error_vec[1] << "\n";
@@ -1078,3 +1081,16 @@ void Tools::PrintErrors(std::ofstream& out, const ProblemConfig& config, const T
     out << ss.str();
     std::cout << ss.str();
 }
+
+void Tools::PrintErrors(std::ostream& out, const TPZVec& error_vec) {
+
+    std::stringstream ss;
+    ss << '\n';
+    ss << "Global estimator = " << error_vec[3] << "\n";
+    ss << "|ufem-urec| = " << error_vec[1] << "\n";
+    ss << "Residual Error L2 = " << error_vec[4] << "\n";
+    ss << "[Unknown exact solution and errors]\n";
+
+    out << ss.str();
+    std::cout << ss.str();
+}
diff --git a/Projects/Tools/Tools.h b/Projects/Tools/Tools.h
index f1de9ea4..aa0b613e 100644
--- a/Projects/Tools/Tools.h
+++ b/Projects/Tools/Tools.h
@@ -14,8 +14,6 @@
 #include "TPZVTKGeoMesh.h"
 #include "ProblemConfig.h"
 
-#include "DarcyFlow/TPZDarcyFlow.h"
-
 #include "pzintel.h"
 
 #include "pzbuildmultiphysicsmesh.h"
@@ -23,7 +21,6 @@
 #include "TPZMultiphysicsInterfaceEl.h"
 #include "TPZHybridizeHDiv.h"
 
-#include "TPZAnalysis.h"
 #include "pzstepsolver.h"
 #include "TPZSSpStructMatrix.h"
 #include "TPZParFrontStructMatrix.h"
@@ -58,7 +55,7 @@ namespace Tools {
 
     TPZCompMesh *CreatePressureMesh(const ProblemConfig &problem);
 
-    TPZMultiphysicsCompMesh *CreateHDivMesh(const ProblemConfig &problem);
+    TPZMultiphysicsCompMesh *CreateMixedMesh(const ProblemConfig &problem);
 
     void UniformRefinement(int nDiv, TPZGeoMesh *gmesh);
 
@@ -82,7 +79,7 @@ namespace Tools {
 
     void SolveMixedProblem(TPZCompMesh *cmesh_HDiv, const ProblemConfig &config);
 
-    TPZCompMesh *CMeshH1(ProblemConfig problem);
+    TPZCompMesh *CMeshH1(const ProblemConfig& problem);
 
     void hAdaptivity(TPZCompMesh *postProcessMesh, TPZGeoMesh *gmeshToRefine, ProblemConfig &config);
 
@@ -97,6 +94,7 @@ namespace Tools {
     TPZGeoMesh* CreateGeoMesh(int nelem, TPZVec& bcids, int dim, bool isOriginCentered, int topologyMode);
 
     void PrintErrors(std::ofstream& out, const ProblemConfig& config, const TPZVec& error_vec);
+    void PrintErrors(std::ostream& out, const TPZVec& error_vec);
 }
 
 #endif
diff --git a/Projects/Tools/ToolsMHM.cpp b/Projects/Tools/ToolsMHM.cpp
index 034824e4..584c182d 100644
--- a/Projects/Tools/ToolsMHM.cpp
+++ b/Projects/Tools/ToolsMHM.cpp
@@ -1,126 +1,24 @@
 //
 //  ToolsMHM.cpp
-//  AdaptivityTest
 //
 //  Created by Denise De Siqueira on 01/11/19.
 //
 
 #include "ToolsMHM.h"
-#include 
 #include "pzvec.h"
 #include "pzstack.h"
-#include "pzfmatrix.h"
 #include "pzreal.h"
 #include "pzgmesh.h"
-#include "pzcmesh.h"
-#include "pzcompel.h"
-#include "pzgeoelside.h"
-#include "tpzgeoelrefpattern.h"
-#include "tpzautopointer.h"
-#include "TPZLinearAnalysis.h"
-
-#include "TPZSSpStructMatrix.h"
-#include "pzstepsolver.h"
-#include "TPZStructMatrixT.h"
-
-#include "tpzarc3d.h"
 #include "tpzgeoblend.h"
-
-
-#include "pzbuildmultiphysicsmesh.h"
-#include "TPZCompMeshTools.h"
-
-#include "TPZVTKGeoMesh.h"
-
 #include "TPZHybridizeHDiv.h"
 #include "TPZGenGrid2D.h"
 
-#include "TPZNullMaterial.h"
-
-#include 
 #include 
 #include 
-#include 
-#include "TPZPersistenceManager.h"
-
-#include "TPZMHMHDivErrorEstimator.h"
-
-TPZCompMesh *CMeshFlux(TPZGeoMesh * gmesh,int pOrder){
-    
-    int impervious_mat = 1;
-    int permeable_mat = 2;
-    int dim = gmesh->Dimension();
-    
-    TPZVec convdir(dim,0.0);
-    
-    TPZCompMesh *cmesh = new TPZCompMesh(gmesh);
-    cmesh->SetDimModel(dim);
-    cmesh->SetDefaultOrder(pOrder); //Default polynomial order of the approximation
-    
-    //Definition of the approximation space:
-    
-    TPZNullMaterial<> *mat_0 = new TPZNullMaterial<>(impervious_mat,dim);
-    TPZNullMaterial<> *mat_1 = new TPZNullMaterial<>(permeable_mat,dim);
-    
-    //  inserting volumetric materials objects
-    cmesh->InsertMaterialObject(mat_0);
-    cmesh->InsertMaterialObject(mat_1);
-    
-    cmesh->SetAllCreateFunctionsHDiv(); //Creating H(div) functions
-    
-    
-    int type_D = 0;
-    int type_N = 1;
-    TPZFMatrix val1(1, 1, 0.);
-    TPZManVector val2(1, 0.);
-    
-    // Insert boundary conditions
-    //Neumann boundary conditions (flux = 0)
-    int right_bc_id = -2;
-    TPZBndCondT * right_bc = mat_0->CreateBC(mat_0, right_bc_id, type_N, val1, val2);
-    cmesh->InsertMaterialObject(right_bc);
-    
-    int left_bc_id = -4;
-    TPZBndCondT * left_bc = mat_0->CreateBC(mat_0, left_bc_id, type_N, val1, val2);
-    cmesh->InsertMaterialObject(left_bc);
-    
-    int bottom_bc_1id = -1;
-    TPZBndCondT* bottom_bc_1 = mat_0->CreateBC(mat_0, bottom_bc_1id, type_N, val1, val2);
-    cmesh->InsertMaterialObject(bottom_bc_1);
-    
-    int top_bc_1id = -3;
-    TPZBndCondT * top_bc_1 = mat_0->CreateBC(mat_0, top_bc_1id, type_N, val1, val2);
-    cmesh->InsertMaterialObject(top_bc_1);
-    
-    
-    //Dirichlet Conditions (p=1 in, p=0 out)
-    int bottom_bc_id = -5;
-    TPZBndCondT * bottom_bc = mat_0->CreateBC(mat_0, bottom_bc_id, type_D, val1, val2);
-    cmesh->InsertMaterialObject(bottom_bc);
-    
-    int top_bc_id = -6;
-    TPZBndCondT * top_bc = mat_0->CreateBC(mat_0, top_bc_id, type_D, val1, val2);
-    cmesh->InsertMaterialObject(top_bc);
-    
-    cmesh->SetName("LaberintoTest");
-    cmesh->AutoBuild();
-    
-#ifdef ERRORESTIMATION_DEBUG
-    std::ofstream file("cmesh_flux.txt");
-    cmesh->Print(file);
-#endif
-    
-    return cmesh;
-    
-}
 
 // compute the coarse indices of the geometric mesh
 void ComputeCoarseIndices(TPZGeoMesh *gmesh, TPZVec &coarseindices)
 {
-    //    {
-    //        std::ofstream out("gmeshref.txt");
-    //        gmesh->Print(out);
-    //    }
     coarseindices.Resize(gmesh->NElements());
     int count = 0;
     for (int64_t el=0; elNElements(); el++) {
@@ -133,225 +31,23 @@ void ComputeCoarseIndices(TPZGeoMesh *gmesh, TPZVec &coarseindices)
     coarseindices.Resize(count);
 }
 
-
-std::map> IdentifyChanel(TPZCompMesh *cmesh){
-    
-    int nelements = cmesh->NElements();
-    TPZGeoElSide first_gelside;
-    for(int iel=0; iel<=nelements; iel++){
-        TPZCompEl *cel = cmesh->Element(iel);
-        if(!cel){continue;}
-        TPZGeoEl *gel = cel->Reference();
-        if(gel->MaterialId() == -5)
-        {
-            TPZGeoElSide gelside(gel);
-            first_gelside = gelside.Neighbour();
-            break;
-        }
-    }
-    if(!first_gelside) DebugStop();
-    int count=0;
-    int count_chain=0;
-    double Flux_Max = 0.0;
-    std::map> chain;
-    bool exit=false;
-    
-    while(exit==false){
-        
-        // if their is a neighbour along side 6 with material id -6, we found the exit
-        TPZGeoElSide exit_test = first_gelside;
-        exit_test.SetSide(6);
-        TPZGeoElSide candidate_exist =exit_test.Neighbour();
-        if(candidate_exist.Element()->MaterialId()==-6){
-            std::cout<<"Cadena encontrada con exito";
-            exit = true;
-            break;
-        }
-        
-        
-        //    while(candidate.Element()->Dimension()!=2){
-        //
-        //            candidate=candidate.Neighbour();
-        //
-        //    }
-        int side_in = first_gelside.Side();
-        // loop over the one dimensional sides
-        for(int ican=4; ican<8; ican++){
-            // if the one-d side is the entry side, do not consider
-            if(side_in == ican ){continue;}
-
-            first_gelside.SetSide(ican);
-            // find a neighbour of dimension 2
-            TPZGeoElSide candidate = first_gelside.Neighbour();
-            while(candidate.Element()->Dimension()!=2){
-                
-                candidate=candidate.Neighbour();
-            }
-            if(candidate.Element() == first_gelside.Element()) DebugStop();
-            
-            //compute the flux values at the center of the 3 neighbours
-            TPZVec qsi(2);
-            qsi[0]=0;
-            qsi[1]=0;
-            int var=1;
-            TPZVec sol;
-            TPZCompEl *cel = candidate.Element()->Reference();
-            cel->Solution(qsi, var,sol);
-            double Flux_can_Mag = sqrt(sol[0]*sol[0] + sol[1]*sol[1]);
-            
-            if(Flux_can_Mag > Flux_Max){
-                Flux_Max =Flux_can_Mag;
-                first_gelside.SetSide(ican);
-                // first is the exit side
-                // second is the entry side
-                chain[count].first= first_gelside;
-                chain[count].second =candidate;
-            }
-        }
-        Flux_Max=0.0;
-        first_gelside =chain[count].second;
-        count++;
-    }
-    
-    //here
-    
-    cmesh->LoadReferences();
-    TPZGeoMesh *gmesh =cmesh->Reference();
-//    for(auto it:gmesh->ElementVec()){
-//        int father_index = it->FatherIndex();
-//        int element = it->Index();
-//        std::cout<<"Element: "<> skelchanel;
-    // the skel channel will filter the elements of chain that
-    // have different father index
-    // create a boundary element of matid 10 on the interfaces between macro
-    // elements
-    int count_skel_chanel=0;
-    int matId_skel_chanel = 10;
-    count =0;
-    for(auto it : chain){
-        TPZGeoEl *first_element = it.second.first.Element();
-        TPZGeoEl *second_element = it.second.second.Element();
-        
-        int first_father_index = first_element->LowestFather()->Index();
-        int second_father_index = second_element->LowestFather()->Index();
-        
-        if(0)
-        {
-            TPZManVector xic1(2),xic2(2),xc1(3),xc2(3);
-            first_element->CenterPoint(8, xic1);
-            first_element->X(xic1, xc1);
-            second_element->CenterPoint(8, xic2);
-            second_element->X(xic2, xc2);
-            std::cout << "fathers " << first_father_index << ' ' << second_father_index << " x1 " << xc1 << " x2 " << xc2 << std::endl;
-        }
-        
-        if(first_father_index!=second_father_index){
-            TPZGeoElBC(it.second.first, matId_skel_chanel);
-            skelchanel[count_skel_chanel] = it.second;
-            count_skel_chanel++;
-        }
-        
-    }
-    std::ofstream out("TestMesh.vtk");
-    TPZVTKGeoMesh::PrintGMeshVTK(gmesh, out, true);
-    
-    return skelchanel;
-    
-}
-
-TPZGeoMesh *CreateCircleGeoMesh() {
-
-    TPZGeoMesh *gmesh = new TPZGeoMesh();
-    gmesh->SetDimension(2);
-
-    TPZVec coord(3, 0.);
-
-    // Inserts node at origin
-    gmesh->NodeVec().AllocateNewElement();
-    gmesh->NodeVec()[0].Initialize(coord, *gmesh);
-
-    // Inserts circumference nodes
-    for (int64_t i = 0; i < 16; i++) {
-        const REAL step = M_PI / 8;
-        coord[0] = cos(i * step);
-        coord[1] = sin(i * step);
-        std::cout<<"{"<NodeVec().AllocateNewElement();
-        gmesh->NodeVec()[newID].Initialize(coord, *gmesh);
-    }
-
-    int matIdTriangle = 1, matIdArc = 2;
-
-    // Inserts triangle elements
-    TPZManVector nodesIdVec(3);
-    for (int64_t i = 0; i < 7; i++) {
-        nodesIdVec[0] = 0;
-        nodesIdVec[1] = 1 + 2 * i;
-        nodesIdVec[2] = 3 + 2 * i;
-        std::cout<<"{"<>(nodesIdVec, matIdTriangle, *gmesh);
-    }
-    nodesIdVec[0] = 0;
-    nodesIdVec[1] = 15;
-    nodesIdVec[2] = 1;
-
-    new TPZGeoElRefPattern>(nodesIdVec, matIdTriangle, *gmesh);
-    // Inserts arc elements
-    for (int64_t i = 0; i < 7; i++) {
-        nodesIdVec[0] = 1 + 2 * i;
-        nodesIdVec[1] = 3 + 2 * i;
-        nodesIdVec[2] = 2 + 2 * i;
-        new TPZGeoElRefPattern(nodesIdVec, matIdArc, *gmesh);
-    }
-
-    nodesIdVec[0] = 15;
-    nodesIdVec[1] = 1;
-    nodesIdVec[2] = 16;
-    //para o broblema do douglas matIdArc=3
-    new TPZGeoElRefPattern(nodesIdVec, 3, *gmesh);
-    //    // Finally, inserts line elements to complete boundary
-    //    nodesIdVec.Resize(2);
-    //    nodesIdVec[0] = 0;
-    //    nodesIdVec[1] = 1;
-    //    new TPZGeoElRefPattern(nodesIdVec, matIdArc, *gmesh);
-    //
-    //    nodesIdVec[0] = 0;
-    //    nodesIdVec[1] = 14;
-    //    new TPZGeoElRefPattern(nodesIdVec, matIdArc, *gmesh);
-
-    gmesh->BuildConnectivity();
-
-    return gmesh;
-}
-
 TPZGeoMesh *CreateLMHMMesh(int nDiv, TPZVec& coarseIndexes) {
+
+    int factor = static_cast(lround(pow(2, nDiv) + 0.5));
     
-    int factor = (int)(pow(2, nDiv) + 0.5);
-    
-    int xElements = factor * 6;
-    int yElements = factor * 4;
+    auto xElements = factor * 6;
+    auto yElements = factor * 4;
     
-    TPZManVector nx(2);
+    TPZManVector nx(2);
     nx[0] = xElements;
     nx[1] = yElements;
     
     TPZManVector x0(3,0.), x1(3,1.);
     x1[2] = 0.;
     
-    
     TPZGenGrid2D gen(nx, x0, x1, 1, 0);
     
-    TPZGeoMesh *gmesh = new TPZGeoMesh();
+    auto *gmesh = new TPZGeoMesh();
     gen.Read(gmesh);
     gen.SetBC(gmesh, 4, -1);
     gen.SetBC(gmesh, 5, -1);
@@ -373,12 +69,12 @@ TPZGeoMesh *CreateLMHMMesh(int nDiv, TPZVec& coarseIndexes) {
             TPZGeoEl *gel = gmesh->ElementVec()[elem];
             if (gel->Dimension() != 2) DebugStop();
             
-            int lineInPattern = elem / nx[0] % 4;
-            int colInPattern = elem % nx[0] % 6;
+            auto lineInPattern = elem / nx[0] % 4;
+            auto colInPattern = elem % nx[0] % 6;
             
             // IDs of elements in the neighbourhood to which a coarse index has been already assigned
-            int leftEl = elem - 1;
-            int bottomEl = elem - nx[0];
+            auto leftEl = elem - 1;
+            auto bottomEl = elem - nx[0];
             
             if (lineInPattern == 0) {
                 if (colInPattern % 2 == 0) {
@@ -416,118 +112,6 @@ TPZGeoMesh *CreateLMHMMesh(int nDiv, TPZVec& coarseIndexes) {
             }
         }
     }
-    
-//    for(int ilinha =0 ; ilinha  cmesh, const TPZVec> &compmeshes,
-                  TPZAnalyticSolution *analytic, const std::string &prefix, TRunConfig config) {
-    //calculo solution
-    bool shouldrenumber = true;
-    TPZLinearAnalysis an(cmesh,shouldrenumber);
-#ifdef PZ_USING_MKL
-    TPZSSpStructMatrix<> strmat(cmesh.operator->());
-    strmat.SetNumThreads(0/*config.n_threads*/);
-#else
-    TPZSkylineStructMatrix strmat(cmesh.operator->());
-    strmat.SetNumThreads(config.n_threads);
-#endif
-
-
-#ifdef ERRORESTIMATION_DEBUG
-    if(0)
-    {
-        std::ofstream file("MeshToSolveProblem.txt");
-        cmesh->Print(file);
-    }
-#endif
-
-
-    an.SetStructuralMatrix(strmat);
-    TPZStepSolver step;
-    step.SetDirect(ELDLt);
-    an.SetSolver(step);
-    std::cout << "Assembling\n";
-    an.Assemble();
-
-    std::cout << "Solving\n";
-    an.Solve();
-    std::cout << "Finished\n";
-    an.LoadSolution(); // compute internal dofs
-
-
-    TPZMFSolutionTransfer transfer;
-    transfer.BuildTransferData(cmesh.operator->());
-//    TPZBuildMultiphysicsMesh::TransferFromMultiPhysics(compmeshes, cmesh);
-    transfer.TransferFromMultiphysics();
-
-    if(0)
-    {
-        std::ofstream out1("mfmesh.txt");
-        cmesh->Print(out1);
-        std::ofstream out2("flux.txt");
-        compmeshes[0]->Print(out2);
-        std::ofstream out3("pressure.txt");
-        compmeshes[1]->Print(out3);
-        std::ofstream out4("transfer.txt");
-        transfer.Print(out4);
-    }
 
-    TPZStack scalnames,vecnames;
-    TPZMaterial *mat = cmesh->FindMaterial(1);
-    if (!mat) {
-        DebugStop();
-    }
-    if (analytic)
-    {
-        an.SetExact(analytic->ExactSolution());
-        scalnames.Push("ExactPressure");
-        vecnames.Push("ExactFlux");
-    }
-    {
-        scalnames.Push("Pressure");
-        //  scalnames.Push("Permeability");
-        vecnames.Push("Flux");
-        //   vecnames.Push("Derivative");
-    }
-
-
-    std::string plotname;
-    {
-        std::stringstream out;
-        out << "MHMHdiv"<<"_kin" <Dimension() , scalnames, vecnames, plotname);
-    an.PostProcess(resolution,cmesh->Dimension() );
-
-
-    if(analytic)
-    {
-        TPZManVector errors(4,0.);
-        an.SetThreadsForError(config.n_threads);
-        an.SetExact(analytic->ExactSolution());
-        an.PostProcessError(errors,false);
-
-        //Erro
-
-        std::ofstream myfile;
-        myfile.open("ArquivosErrosMHM.txt", std::ios::app);
-        myfile << "\n\n Error for MHM formulation " ;
-        myfile << "\n-------------------------------------------------- \n";
-        myfile << "Ndiv = " << config.numHDivisions << " Order Internal= " << config.pOrderInternal <<" Order Skeleton= " << config.pOrderSkeleton <<"\n";
-        myfile << "DOF Total = " << cmesh->NEquations() << "\n";
-        myfile << "Energy norm = " << errors[0] << "\n";//norma energia
-        myfile << "error norm L2 = " << errors[1] << "\n";//norma L2
-        myfile << "Semi norm H1 = " << errors[2] << "\n";//norma L2
-        myfile.close();
-
-    }
-
-}
+    return gmesh;
+}
\ No newline at end of file
diff --git a/Projects/Tools/ToolsMHM.h b/Projects/Tools/ToolsMHM.h
index 68f2abf3..b27f34c0 100644
--- a/Projects/Tools/ToolsMHM.h
+++ b/Projects/Tools/ToolsMHM.h
@@ -1,162 +1,18 @@
 //
 //  ToolsMHM.hpp
-//  AdaptivityTest
 //
 //  Created by Denise De Siqueira on 01/11/19.
 //
 
-#ifndef ToolsMHM_hpp
-#define ToolsMHM_hpp
+#ifndef TOOLSMHM_H
+#define TOOLSMHM_H
 
-#include 
-
-#endif /* ToolsMHM_hpp */
-
-
-#ifdef HAVE_CONFIG_H
-#include 
-#endif
-
-#include 
-#include 
 #include "pzvec.h"
-#include "pzstack.h"
-#include "pzfmatrix.h"
-#include "pzfstrmatrix.h"
-#include "pzlog.h"
-#include "pzreal.h"
 #include "pzgmesh.h"
-#include "pzcmesh.h"
-#include "pzcompel.h"
-#include "TPZInterfaceEl.h"
-#include "pzgeoelside.h"
-#include "TPZGeoLinear.h"
-#include "pzgeopoint.h"
-
-#include "TPZRefPattern.h"
-#include "tpzgeoelrefpattern.h"
-#include "tpzautopointer.h"
-#include "TPZAnalysis.h"
-
-#include "TPZSSpStructMatrix.h"
-#include "pzstepsolver.h"
-#include "TPZStructMatrixT.h"
-#include "pzfstrmatrix.h"
-#include "TPZFrontNonSym.h"
-#include "TPZFrontSym.h"
-#include "TPZSpStructMatrix.h"
-#include "pzbstrmatrix.h"
-
-#include "tpzarc3d.h"
-#include "tpzgeoblend.h"
-#include "TPZGeoLinear.h"
-
-#include "DarcyFlow/TPZDarcyFlow.h"
-#include "TPZMatLaplacianHybrid.h"
-#include "TPZLagrangeMultiplierCS.h"
-
-#include "pzbuildmultiphysicsmesh.h"
-#include "pzelementgroup.h"
-#include "TPZCompMeshTools.h"
-#include "pzcondensedcompel.h"
-#include "pzfunction.h"
-#include "pzgraphmesh.h"
-#include "pzfmatrix.h"
-
-#include "pzlog.h"
-
-#include "TPZVTKGeoMesh.h"
-#include "pzvisualmatrix.h"
-
-//#include "pzgengrid.h"
-
-#include "TPZGenGrid2D.h"
-#include "TPZExtendGridDimension.h"
-#include "pzcheckgeom.h"
-
-#include "TPZMHMeshControl.h"
-#include "TPZMHMixedMeshControl.h"
-#include "TPZMHMixedMeshChannelControl.h"
-
-#include "TPZMHMixedHybridMeshControl.h"
-#include "TPZHybridizeHDiv.h"
-//#include "ConfigCasesMaze.h"
-#include 
-#include 
-#include 
-#include 
-
-#include "TPZPersistenceManager.h"
-
-#include "TPZMHMHDivErrorEstimator.h"
-#include "Tools.h"
-#include "TPZGenGrid2D.h"
-
-struct TRunConfig
-{
-    int nelxcoarse = -1;
-    int nelycoarse = -1;
-    int numHDivisions = 0;
-    int pOrderInternal = 1;
-    int hdivmaismais = 1;
-    int numDivSkeleton = 0;
-    int pOrderSkeleton = 1;
-    int Hybridize = 0;
-    int Condensed = 1;
-    int LagrangeMult = 0;
-    int newline = 0;
-    int n_threads = 0;
-    bool MHM_HDiv_Elast = false;
-
-    /// number of equations when not condensing anything
-    int64_t fGlobalSystemSize = -1;
-    /// number of equations considering local condensation
-    int64_t fGlobalSystemWithLocalCondensationSize = -1;
-    /// number of equations of the global system
-    int64_t fNumeq = -1;
-
-    REAL fDeltaT = 1.;
-
-    /// number of timesteps
-    int64_t nTimeSteps = 10;
-
-    std::ostream &InlinePrint(std::ostream &out)
-    {
-//        out << "nelxCoarse " << nelxcoarse << " nelyCoarse " << nelycoarse << " numHDiv " << numHDivisions << " porderInternal " << pOrderInternal << " numDivSkeleton " << numDivSkeleton
-//        << " porderSkeleton " << pOrderSkeleton << " Hybridize " << Hybridize << " Condensed " << Condensed << " LagrangeMult " << LagrangeMult
-//        << " sysnocondense " << fGlobalSystemSize << " syslocalcondense " << fGlobalSystemWithLocalCondensationSize << " neq " << fNumeq;
-        // : n_dx n_dy k_skeleton m_div k_subelement
-        out << "n_dx " << nelxcoarse << " n_dy " << nelycoarse  << " k_skeleton " << pOrderSkeleton << " m_div " << numHDivisions << " k_subelement " << pOrderInternal << " upscaling_dof " << fNumeq << " total_dof " << fGlobalSystemSize;
-        return out;
-    }
-    std::ostream &MathematicaInlinePrint(std::ostream &out)
-    {
-        out << "nelxCoarse, " << nelxcoarse << ", nelyCoarse, " << nelycoarse << " ,numHDiv, " << numHDivisions << " ,porderInternal, " << pOrderInternal << " ,numDivSkeleton, " << numDivSkeleton
-            << " ,porderSkeleton, " << pOrderSkeleton << " ,Hybridize, " << Hybridize << " ,Condensed, " << Condensed << " ,LagrangeMult, " << LagrangeMult
-            << " ,sysnocondense, " << fGlobalSystemSize << " ,syslocalcondense, " << fGlobalSystemWithLocalCondensationSize << " ,neq, " << fNumeq;
-        return out;
-    }
-
-    std::ostream &ConfigPrint(std::ostream &out)
-    {
-        out << nelxcoarse << "x" << nelycoarse << "_HSkel" << numDivSkeleton << "_pSkel" << pOrderSkeleton << "_HDiv" << numHDivisions << "_pInt" << pOrderInternal;
-        return out;
-    }
-};
-// Creating the computational flux mesh
-TPZCompMesh *CMeshFlux(TPZGeoMesh * gmesh,int pOrder);
-
-// Create a geometric mesh with the given parameters, nx and ny are the coarse elements number. The total number of elements are defined by the image read.
-TPZGeoMesh *GenerateGeoMesh(std::string name, int nx, int ny);
-
-// Create a geoElSide map to be open because it has flux.
-// cmesh is the flux mesh obtained from a Hdiv previus simulation.
-std::map> IdentifyChanel (TPZCompMesh *cmesh);
 
 // Compute the geometric mesh coarse indices
 void ComputeCoarseIndices(TPZGeoMesh *gmesh, TPZVec &coarseindices);
 
-TPZGeoMesh *CreateLMHMMesh(int nDiv, TPZVec& coarseIndexes);
+TPZGeoMesh *CreateLMHMMesh(int nDiv, TPZVec &coarseIndexes);
 
-/// Solve the problem composed of a multiphysics mesh composed of compmeshes - applies to MHM and MHM-H(div)
-void SolveProblem(TPZAutoPointer cmesh, const TPZVec >& compmeshes, TPZAnalyticSolution *analytic, const std::string& prefix, TRunConfig config);
+#endif // TOOLSMHM_H