Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21761) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21762) @@ -39,16 +39,16 @@ AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ - bool ismismip = true; - if(ismismip){//itapopo - TPZFileStream fstr; - std::stringstream ss; + //bool ismismip = false; + //if(ismismip){//itapopo + // TPZFileStream fstr; + // std::stringstream ss; - ss << this->levelmax; - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; + // ss << this->levelmax; + // std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; - fstr.OpenWrite(AMRfile.c_str()); - int withclassid = 1; - this->Write(fstr,withclassid); - } + // fstr.OpenWrite(AMRfile.c_str()); + // int withclassid = 1; + // this->Write(fstr,withclassid); + //} this->CleanUp(); gRefDBase.clear(); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21761) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21762) @@ -1621,5 +1621,14 @@ name==VelEnum || name==VxPicardEnum || - name==VyPicardEnum + name==VyPicardEnum || + name==DeviatoricStressxxEnum || + name==DeviatoricStressyyEnum || + name==DeviatoricStressxyEnum || + name==DeviatoricStressxzEnum || + name==DeviatoricStressyzEnum || + name==DeviatoricStresszzEnum || + name==DeviatoricStresseffectiveEnum || + name==VxAverageEnum || + name==VyAverageEnum ) { Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21761) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21762) @@ -2956,11 +2956,22 @@ this->SetRefPatterns(); if(my_rank==0){ - bool ismisomip = true; + bool ismisomip = false; if(ismisomip){//itapopo TPZFileStream fstr; std::stringstream ss; - - ss << levelmax; - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; + int frictionlaw; + + this->parameters->FindParam(&frictionlaw,FrictionLawEnum); + + ss << levelmax; + if(frictionlaw==1){ + ss << "_viscous/amr.txt"; + }else if(frictionlaw==7){ + ss << "_tsai/amr.txt"; + }else{ + _error_("friction law not supported here."); + } + + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); fstr.OpenRead(AMRfile.c_str()); @@ -3135,5 +3146,5 @@ /*Insert MISMIP+ bed topography*/ - if(true) this->BedrockFromMismipPlus(); + if(false) this->BedrockFromMismipPlus(); /*Adjust base, thickness and mask grounded ice leve set*/ @@ -3278,4 +3289,26 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); + + //itapopo + if(IssmComm::GetRank()==0){ + TPZFileStream fstr; + std::stringstream ss; + int frictionlaw,levelmax; + this->parameters->FindParam(&frictionlaw,FrictionLawEnum); + this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); + ss << levelmax; + if(frictionlaw==1){ + ss << "_viscous/amr.txt"; + }else if(frictionlaw==7){ + ss << "_tsai/amr.txt"; + }else{ + _error_("friction law not supported here."); + } + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); + fstr.OpenWrite(AMRfile.c_str()); + int withclassid = 1; + this->amr->Write(fstr,withclassid); + } + //itapopo /*Cleanup*/ @@ -3496,19 +3529,26 @@ /*elements is in Matlab indexing*/ - - int my_rank = IssmComm::GetRank(); - int numberofsegments = -1; - IssmDouble* vx = NULL; //itapopo this is not being used - IssmDouble* vy = NULL; //itapopo this is not being used - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elementslist = NULL; - int* segments = NULL; - IssmDouble* masklevelset = NULL; - const int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements + int my_rank = IssmComm::GetRank(); + int numberofsegments = -1; + IssmDouble* vx = NULL; //itapopo this is not being used + IssmDouble* vy = NULL; //itapopo this is not being used + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elementslist = NULL; + int* segments = NULL; + IssmDouble* masklevelset = NULL; + IssmDouble* pelementerror = NULL; + const int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements /*Solutions which will be used to refine the elements*/ this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse + + //Compute the ZZ error estimator per element + this->ZZErrorEstimator(&pelementerror); + + _printf0_("P Element error\n"); + for(int i=0;ielements->NumberOfElements();i++) _printf0_(""<(segments); xDelete(masklevelset); + xDelete(pelementerror); } @@ -3548,5 +3589,4 @@ int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements - int numberofelements = this->elements->NumberOfElements(); int numberofvertices = this->vertices->NumberOfVertices(); @@ -4006,3 +4046,165 @@ } /*}}}*/ +void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ + + this->DeviatoricStressx();//itapopo + + int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements + int numberofvertices = this->vertices->NumberOfVertices(); + + IssmDouble weight = 0.; + IssmDouble* tauxx = NULL; + IssmDouble* tauyy = NULL; + IssmDouble* tauxy = NULL; + IssmDouble* totalweight = NULL; + IssmDouble* deviatoricstressxx = xNew(elementswidth); + IssmDouble* deviatoricstressyy = xNew(elementswidth); + IssmDouble* deviatoricstressxy = xNew(elementswidth); + int* elem_vertices = xNew(elementswidth); + Vector* vectauxx = new Vector(numberofvertices); + Vector* vectauyy = new Vector(numberofvertices); + Vector* vectauxy = new Vector(numberofvertices); + Vector* vectotalweight = new Vector(numberofvertices); + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element->GetInputListOnVertices(deviatoricstressxx,DeviatoricStressxxEnum); + element->GetInputListOnVertices(deviatoricstressyy,DeviatoricStressyyEnum); + element->GetInputListOnVertices(deviatoricstressxy,DeviatoricStressxyEnum); + element->GetVerticesSidList(elem_vertices); + + /*weight to calculate the smoothed deviatoric stress*/ + Tria* triaelement = xDynamicCast(element); + weight = triaelement->GetArea();//the tria area is a choice for the weight + + /*taux xx*/ + vectauxx->SetValue(elem_vertices[0],weight*deviatoricstressxx[0],ADD_VAL); + vectauxx->SetValue(elem_vertices[1],weight*deviatoricstressxx[1],ADD_VAL); + vectauxx->SetValue(elem_vertices[2],weight*deviatoricstressxx[2],ADD_VAL); + /*tau yy*/ + vectauyy->SetValue(elem_vertices[0],weight*deviatoricstressyy[0],ADD_VAL); + vectauyy->SetValue(elem_vertices[1],weight*deviatoricstressyy[1],ADD_VAL); + vectauyy->SetValue(elem_vertices[2],weight*deviatoricstressyy[2],ADD_VAL); + /*tau xy*/ + vectauxy->SetValue(elem_vertices[0],weight*deviatoricstressxy[0],ADD_VAL); + vectauxy->SetValue(elem_vertices[1],weight*deviatoricstressxy[1],ADD_VAL); + vectauxy->SetValue(elem_vertices[2],weight*deviatoricstressxy[2],ADD_VAL); + /*total weight*/ + vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL); + vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL); + vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL); + } + + /*Assemble*/ + vectauxx->Assemble(); + vectauyy->Assemble(); + vectauxy->Assemble(); + vectotalweight->Assemble(); + + /*Serialize*/ + tauxx = vectauxx->ToMPISerial(); + tauyy = vectauyy->ToMPISerial(); + tauxy = vectauxy->ToMPISerial(); + totalweight = vectotalweight->ToMPISerial(); + + /*Divide for the total weight*/ + for(int i=0;i0); + tauxx[i] = tauxx[i]/totalweight[i]; + tauyy[i] = tauyy[i]/totalweight[i]; + tauxy[i] = tauxy[i]/totalweight[i]; + } + + /*Set output*/ + (*ptauxx) = tauxx; + (*ptauyy) = tauyy; + (*ptauxy) = tauxy; + + /*Cleanup*/ + delete vectauxx; + delete vectauyy; + delete vectauxy; + delete vectotalweight; + xDelete(deviatoricstressxx); + xDelete(deviatoricstressyy); + xDelete(deviatoricstressxy); + xDelete(totalweight); + xDelete(elem_vertices); + +} +/*}}}*/ +void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ + + /*Compute the Zienkiewicz and Zhu (ZZ) error estimator. + * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ + + IssmDouble Jdet,error,ftxx,ftyy,ftxy; + int sid; + int numnodes = this->GetElementsWidth();//just 2D mesh, tria elements, P1 + int numberofelements = this->elements->NumberOfElements(); + IssmDouble* xyz_list = NULL; + IssmDouble* smoothedtauxx = NULL; + IssmDouble* smoothedtauyy = NULL; + IssmDouble* smoothedtauxy = NULL; + IssmDouble* tauxx = xNew(numnodes); + IssmDouble* tauyy = xNew(numnodes); + IssmDouble* tauxy = xNew(numnodes); + IssmDouble* basis = xNew(numnodes); + int* elem_vertices = xNew(numnodes); + Vector* velementerror= new Vector(numberofelements); + + /*Get smoothed deviatoric stress tensor*/ + this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); + + /*Integrate the error over elements*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element->GetInputListOnVertices(tauxx,DeviatoricStressxxEnum); + element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum); + element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); + element->GetVerticesSidList(elem_vertices); + + /*Integrate*/ + element->GetVerticesCoordinates(&xyz_list); + Gauss* gauss=element->NewGauss(2); + error=0.; + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + ftxx=0;ftyy=0;ftxy=0; + for(int n=0;nweight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) ); + } + /*Set the error in the global vector*/ + sid=element->Sid(); + velementerror->SetValue(sid,error,INS_VAL); + /*Cleanup intermediaries*/ + xDelete(xyz_list); + delete gauss; + } + + /*Assemble*/ + velementerror->Assemble(); + + /*Serialize and set output*/ + (*pelementerror)=velementerror->ToMPISerial(); + + /*Cleanup*/ + xDelete(smoothedtauxx); + xDelete(smoothedtauyy); + xDelete(smoothedtauxy); + xDelete(tauxx); + xDelete(tauyy); + xDelete(tauxy); + xDelete(basis); + xDelete(elem_vertices); + delete velementerror; + +} +/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21761) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21762) @@ -166,4 +166,6 @@ void WriteMeshInResults(void); void SetRefPatterns(void); + void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY + void ZZErrorEstimator(IssmDouble** pelementerror); #endif };