![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21762
- Timestamp:
- 06/12/17 04:59:10 (8 years ago)
- Location:
- issm/trunk-jpl/src/c/classes
- Files:
-
- 4 edited
-
AdaptiveMeshRefinement.cpp (modified) (1 diff)
-
Elements/Element.cpp (modified) (1 diff)
-
FemModel.cpp (modified) (7 diffs)
-
FemModel.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
r21672 r21762 39 39 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ 40 40 41 bool ismismip = true;42 if(ismismip){//itapopo43 TPZFileStream fstr;44 std::stringstream ss;41 //bool ismismip = false; 42 //if(ismismip){//itapopo 43 // TPZFileStream fstr; 44 // std::stringstream ss; 45 45 46 ss << this->levelmax;47 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt";46 // ss << this->levelmax; 47 // std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; 48 48 49 fstr.OpenWrite(AMRfile.c_str());50 int withclassid = 1;51 this->Write(fstr,withclassid);52 }49 // fstr.OpenWrite(AMRfile.c_str()); 50 // int withclassid = 1; 51 // this->Write(fstr,withclassid); 52 //} 53 53 this->CleanUp(); 54 54 gRefDBase.clear(); -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r21709 r21762 1621 1621 name==VelEnum || 1622 1622 name==VxPicardEnum || 1623 name==VyPicardEnum 1623 name==VyPicardEnum || 1624 name==DeviatoricStressxxEnum || 1625 name==DeviatoricStressyyEnum || 1626 name==DeviatoricStressxyEnum || 1627 name==DeviatoricStressxzEnum || 1628 name==DeviatoricStressyzEnum || 1629 name==DeviatoricStresszzEnum || 1630 name==DeviatoricStresseffectiveEnum || 1631 name==VxAverageEnum || 1632 name==VyAverageEnum 1624 1633 1625 1634 ) { -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r21706 r21762 2956 2956 this->SetRefPatterns(); 2957 2957 if(my_rank==0){ 2958 bool ismisomip = true;2958 bool ismisomip = false; 2959 2959 if(ismisomip){//itapopo 2960 2960 TPZFileStream fstr; 2961 2961 std::stringstream ss; 2962 2963 ss << levelmax; 2964 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt"; 2962 int frictionlaw; 2963 2964 this->parameters->FindParam(&frictionlaw,FrictionLawEnum); 2965 2966 ss << levelmax; 2967 if(frictionlaw==1){ 2968 ss << "_viscous/amr.txt"; 2969 }else if(frictionlaw==7){ 2970 ss << "_tsai/amr.txt"; 2971 }else{ 2972 _error_("friction law not supported here."); 2973 } 2974 2975 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); 2965 2976 fstr.OpenRead(AMRfile.c_str()); 2966 2977 … … 3135 3146 3136 3147 /*Insert MISMIP+ bed topography*/ 3137 if( true) this->BedrockFromMismipPlus();3148 if(false) this->BedrockFromMismipPlus(); 3138 3149 3139 3150 /*Adjust base, thickness and mask grounded ice leve set*/ … … 3278 3289 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum, 3279 3290 y,numberofvertices,1,step,time)); 3291 3292 //itapopo 3293 if(IssmComm::GetRank()==0){ 3294 TPZFileStream fstr; 3295 std::stringstream ss; 3296 int frictionlaw,levelmax; 3297 this->parameters->FindParam(&frictionlaw,FrictionLawEnum); 3298 this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); 3299 ss << levelmax; 3300 if(frictionlaw==1){ 3301 ss << "_viscous/amr.txt"; 3302 }else if(frictionlaw==7){ 3303 ss << "_tsai/amr.txt"; 3304 }else{ 3305 _error_("friction law not supported here."); 3306 } 3307 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); 3308 fstr.OpenWrite(AMRfile.c_str()); 3309 int withclassid = 1; 3310 this->amr->Write(fstr,withclassid); 3311 } 3312 //itapopo 3280 3313 3281 3314 /*Cleanup*/ … … 3496 3529 3497 3530 /*elements is in Matlab indexing*/ 3498 3499 int my_rank = IssmComm::GetRank();3500 int numberofsegments = -1;3501 IssmDouble* v x= NULL; //itapopo this is not being used3502 IssmDouble* vy = NULL; //itapopo this is not being used3503 IssmDouble* x= NULL;3504 IssmDouble* y= NULL;3505 IssmDouble* z= NULL;3506 int* elementslist= NULL;3507 int* segments= NULL;3508 IssmDouble* masklevelset= NULL;3509 const int elementswidth 3531 int my_rank = IssmComm::GetRank(); 3532 int numberofsegments = -1; 3533 IssmDouble* vx = NULL; //itapopo this is not being used 3534 IssmDouble* vy = NULL; //itapopo this is not being used 3535 IssmDouble* x = NULL; 3536 IssmDouble* y = NULL; 3537 IssmDouble* z = NULL; 3538 int* elementslist = NULL; 3539 int* segments = NULL; 3540 IssmDouble* masklevelset = NULL; 3541 IssmDouble* pelementerror = NULL; 3542 const int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements 3510 3543 3511 3544 /*Solutions which will be used to refine the elements*/ 3512 3545 this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse 3546 3547 //Compute the ZZ error estimator per element 3548 this->ZZErrorEstimator(&pelementerror); 3549 3550 _printf0_("P Element error\n"); 3551 for(int i=0;i<this->elements->NumberOfElements();i++) _printf0_(""<<pelementerror[i]<< "\n"); 3552 _printf0_("\n"); 3513 3553 3514 3554 if(my_rank==0){ … … 3542 3582 if(segments) xDelete<int>(segments); 3543 3583 xDelete<IssmDouble>(masklevelset); 3584 xDelete<IssmDouble>(pelementerror); 3544 3585 3545 3586 } … … 3548 3589 3549 3590 int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements 3550 int numberofelements = this->elements->NumberOfElements();3551 3591 int numberofvertices = this->vertices->NumberOfVertices(); 3552 3592 … … 4006 4046 } 4007 4047 /*}}}*/ 4048 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ 4049 4050 this->DeviatoricStressx();//itapopo 4051 4052 int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements 4053 int numberofvertices = this->vertices->NumberOfVertices(); 4054 4055 IssmDouble weight = 0.; 4056 IssmDouble* tauxx = NULL; 4057 IssmDouble* tauyy = NULL; 4058 IssmDouble* tauxy = NULL; 4059 IssmDouble* totalweight = NULL; 4060 IssmDouble* deviatoricstressxx = xNew<IssmDouble>(elementswidth); 4061 IssmDouble* deviatoricstressyy = xNew<IssmDouble>(elementswidth); 4062 IssmDouble* deviatoricstressxy = xNew<IssmDouble>(elementswidth); 4063 int* elem_vertices = xNew<int>(elementswidth); 4064 Vector<IssmDouble>* vectauxx = new Vector<IssmDouble>(numberofvertices); 4065 Vector<IssmDouble>* vectauyy = new Vector<IssmDouble>(numberofvertices); 4066 Vector<IssmDouble>* vectauxy = new Vector<IssmDouble>(numberofvertices); 4067 Vector<IssmDouble>* vectotalweight = new Vector<IssmDouble>(numberofvertices); 4068 4069 for(int i=0;i<this->elements->Size();i++){ 4070 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 4071 element->GetInputListOnVertices(deviatoricstressxx,DeviatoricStressxxEnum); 4072 element->GetInputListOnVertices(deviatoricstressyy,DeviatoricStressyyEnum); 4073 element->GetInputListOnVertices(deviatoricstressxy,DeviatoricStressxyEnum); 4074 element->GetVerticesSidList(elem_vertices); 4075 4076 /*weight to calculate the smoothed deviatoric stress*/ 4077 Tria* triaelement = xDynamicCast<Tria*>(element); 4078 weight = triaelement->GetArea();//the tria area is a choice for the weight 4079 4080 /*taux xx*/ 4081 vectauxx->SetValue(elem_vertices[0],weight*deviatoricstressxx[0],ADD_VAL); 4082 vectauxx->SetValue(elem_vertices[1],weight*deviatoricstressxx[1],ADD_VAL); 4083 vectauxx->SetValue(elem_vertices[2],weight*deviatoricstressxx[2],ADD_VAL); 4084 /*tau yy*/ 4085 vectauyy->SetValue(elem_vertices[0],weight*deviatoricstressyy[0],ADD_VAL); 4086 vectauyy->SetValue(elem_vertices[1],weight*deviatoricstressyy[1],ADD_VAL); 4087 vectauyy->SetValue(elem_vertices[2],weight*deviatoricstressyy[2],ADD_VAL); 4088 /*tau xy*/ 4089 vectauxy->SetValue(elem_vertices[0],weight*deviatoricstressxy[0],ADD_VAL); 4090 vectauxy->SetValue(elem_vertices[1],weight*deviatoricstressxy[1],ADD_VAL); 4091 vectauxy->SetValue(elem_vertices[2],weight*deviatoricstressxy[2],ADD_VAL); 4092 /*total weight*/ 4093 vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL); 4094 vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL); 4095 vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL); 4096 } 4097 4098 /*Assemble*/ 4099 vectauxx->Assemble(); 4100 vectauyy->Assemble(); 4101 vectauxy->Assemble(); 4102 vectotalweight->Assemble(); 4103 4104 /*Serialize*/ 4105 tauxx = vectauxx->ToMPISerial(); 4106 tauyy = vectauyy->ToMPISerial(); 4107 tauxy = vectauxy->ToMPISerial(); 4108 totalweight = vectotalweight->ToMPISerial(); 4109 4110 /*Divide for the total weight*/ 4111 for(int i=0;i<numberofvertices;i++){ 4112 _assert_(totalweight[i]>0); 4113 tauxx[i] = tauxx[i]/totalweight[i]; 4114 tauyy[i] = tauyy[i]/totalweight[i]; 4115 tauxy[i] = tauxy[i]/totalweight[i]; 4116 } 4117 4118 /*Set output*/ 4119 (*ptauxx) = tauxx; 4120 (*ptauyy) = tauyy; 4121 (*ptauxy) = tauxy; 4122 4123 /*Cleanup*/ 4124 delete vectauxx; 4125 delete vectauyy; 4126 delete vectauxy; 4127 delete vectotalweight; 4128 xDelete<IssmDouble>(deviatoricstressxx); 4129 xDelete<IssmDouble>(deviatoricstressyy); 4130 xDelete<IssmDouble>(deviatoricstressxy); 4131 xDelete<IssmDouble>(totalweight); 4132 xDelete<int>(elem_vertices); 4133 4134 } 4135 /*}}}*/ 4136 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ 4137 4138 /*Compute the Zienkiewicz and Zhu (ZZ) error estimator. 4139 * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ 4140 4141 IssmDouble Jdet,error,ftxx,ftyy,ftxy; 4142 int sid; 4143 int numnodes = this->GetElementsWidth();//just 2D mesh, tria elements, P1 4144 int numberofelements = this->elements->NumberOfElements(); 4145 IssmDouble* xyz_list = NULL; 4146 IssmDouble* smoothedtauxx = NULL; 4147 IssmDouble* smoothedtauyy = NULL; 4148 IssmDouble* smoothedtauxy = NULL; 4149 IssmDouble* tauxx = xNew<IssmDouble>(numnodes); 4150 IssmDouble* tauyy = xNew<IssmDouble>(numnodes); 4151 IssmDouble* tauxy = xNew<IssmDouble>(numnodes); 4152 IssmDouble* basis = xNew<IssmDouble>(numnodes); 4153 int* elem_vertices = xNew<int>(numnodes); 4154 Vector<IssmDouble>* velementerror= new Vector<IssmDouble>(numberofelements); 4155 4156 /*Get smoothed deviatoric stress tensor*/ 4157 this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); 4158 4159 /*Integrate the error over elements*/ 4160 for(int i=0;i<this->elements->Size();i++){ 4161 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 4162 element->GetInputListOnVertices(tauxx,DeviatoricStressxxEnum); 4163 element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum); 4164 element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); 4165 element->GetVerticesSidList(elem_vertices); 4166 4167 /*Integrate*/ 4168 element->GetVerticesCoordinates(&xyz_list); 4169 Gauss* gauss=element->NewGauss(2); 4170 error=0.; 4171 for(int ig=gauss->begin();ig<gauss->end();ig++){ 4172 gauss->GaussPoint(ig); 4173 element->JacobianDeterminant(&Jdet,xyz_list,gauss); 4174 element->NodalFunctions(basis,gauss); 4175 ftxx=0;ftyy=0;ftxy=0; 4176 for(int n=0;n<numnodes;n++) { 4177 ftxx+=(tauxx[n]-smoothedtauxx[elem_vertices[n]])*basis[n]; 4178 ftyy+=(tauyy[n]-smoothedtauyy[elem_vertices[n]])*basis[n]; 4179 ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n]; 4180 } 4181 error+=Jdet*gauss->weight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) ); 4182 } 4183 /*Set the error in the global vector*/ 4184 sid=element->Sid(); 4185 velementerror->SetValue(sid,error,INS_VAL); 4186 /*Cleanup intermediaries*/ 4187 xDelete<IssmDouble>(xyz_list); 4188 delete gauss; 4189 } 4190 4191 /*Assemble*/ 4192 velementerror->Assemble(); 4193 4194 /*Serialize and set output*/ 4195 (*pelementerror)=velementerror->ToMPISerial(); 4196 4197 /*Cleanup*/ 4198 xDelete<IssmDouble>(smoothedtauxx); 4199 xDelete<IssmDouble>(smoothedtauyy); 4200 xDelete<IssmDouble>(smoothedtauxy); 4201 xDelete<IssmDouble>(tauxx); 4202 xDelete<IssmDouble>(tauyy); 4203 xDelete<IssmDouble>(tauxy); 4204 xDelete<IssmDouble>(basis); 4205 xDelete<int>(elem_vertices); 4206 delete velementerror; 4207 4208 } 4209 /*}}}*/ 4008 4210 #endif -
issm/trunk-jpl/src/c/classes/FemModel.h
r21706 r21762 166 166 void WriteMeshInResults(void); 167 167 void SetRefPatterns(void); 168 void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY 169 void ZZErrorEstimator(IssmDouble** pelementerror); 168 170 #endif 169 171 };
Note:
See TracChangeset
for help on using the changeset viewer.
