Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 20116) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 20117) @@ -257,4 +257,5 @@ ./cores/damage_core.cpp\ ./cores/levelsetfunctionslope_core.cpp\ + ./cores/levelset_core.cpp\ ./modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp\ ./classes/Loads/Riftfront.cpp\ Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 20116) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 20117) @@ -59,5 +59,6 @@ bool save_results; int extvar_enum; - femmodel->parameters->FindParam(&extvar_enum, ExtrapolationVariableEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&extvar_enum, ExtrapolationVariableEnum); /*activate formulation: */ @@ -67,9 +68,7 @@ solutionsequence_linear(femmodel); - save_results=true; if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - int outputs[2] = {VxEnum,VyEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + femmodel->RequestedOutputsx(&femmodel->results,&extvar_enum,1); } }/*}}}*/ @@ -84,30 +83,36 @@ ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/ - if(!element->IsOnBase()) return NULL; - Element* basalelement = element->SpawnBasalElement(); - /*Intermediaries */ - int domaintype,dim; - int i,row,col,stabilization; - bool extrapolatebydiffusion = true; + int dim, domaintype, extrapolationcase; + int i,row,col,stabilization; + bool extrapolatebydiffusion; IssmDouble Jdet,D_scalar,h; IssmDouble norm_dlsf; IssmDouble hx,hy,hz,kappa; - IssmDouble* xyz_list = NULL; - - /*Get problem dimension*/ - basalelement->FindParam(&domaintype,DomainTypeEnum); + IssmDouble* xyz_list = NULL; + Element* workelement=NULL; + + /*Get problem case*/ + extrapolationcase=GetExtrapolationCase(element); + switch(extrapolationcase){ + case 0: + if(!element->IsOnBase()) return NULL; + workelement = element->SpawnBasalElement(); + break; + case 1: case 2: case 3: workelement=element; break; + } + /* get extrapolation dimension */ + workelement->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DverticalEnum: dim = 1; break; - case Domain2DhorizontalEnum: dim = 2; break; - case Domain3DEnum: dim = 2; break; - default: _error_("mesh "<GetNumberOfNodes(); + int numnodes = workelement->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = basalelement->NewElementMatrix(); + ElementMatrix* Ke = workelement->NewElementMatrix(); IssmDouble* B = xNew(dim*numnodes); IssmDouble* Bprime = xNew(dim*numnodes); @@ -117,25 +122,25 @@ /*Retrieve all inputs and parameters*/ - Input* lsf_slopex_input=basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - Input* lsf_slopey_input=basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); - basalelement->GetVerticesCoordinates(&xyz_list); - h = basalelement->CharacteristicLength(); + Input* lsf_slopex_input=workelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + Input* lsf_slopey_input=workelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + workelement->GetVerticesCoordinates(&xyz_list); /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss=workelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){/*{{{*/ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - GetB(B,basalelement,xyz_list,gauss); - GetBprime(Bprime,basalelement,xyz_list,gauss); + workelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + GetB(B,workelement,xyz_list,gauss,dim); + GetBprime(Bprime,workelement,xyz_list,gauss,dim); D_scalar=gauss->weight*Jdet; + extrapolatebydiffusion=true; if(extrapolatebydiffusion){ /* diffuse values outward */ for(row=0;rowGetInputValue(&dlsf[0],gauss); - lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + if(dim>1) + lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + if(dim>2) + dlsf[2]=0.; norm_dlsf=0.; for(i=0;ivalues[0],1); - /* Stabilization *//*{{{*/ + /* stabilization *//*{{{*/ + /* do not use streamline upwinding for extrapolation: it yields oscillating results due to diffusion along normal direction, but none across */ stabilization=1; if (stabilization==0){/* no stabilization, do nothing*/} else if(stabilization==1){ /* Artificial Diffusion */ - basalelement->ElementSizes(&hx,&hy,&hz); - h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2)); + workelement->ElementSizes(&hx,&hy,&hz); + h=sqrt(pow(hx*normal[0],2) + pow(hy*normal[1],2)); kappa=h/2.+1.e-14; for(row=0;rowvalues[0],1); } - else if(stabilization==2){ - /*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */ - for(row=0;rowvalues[0],1); - }/*}}}*/ + /*}}}*/ } }/*}}}*/ @@ -213,5 +212,5 @@ xDelete(normal); delete gauss; - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;}; return Ke; @@ -219,23 +218,29 @@ ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/ - if(!element->IsOnBase()) return NULL; - Element* basalelement = element->SpawnBasalElement(); - /*Intermediaries */ - int i, domaintype; - + Element* workelement=NULL; + + /*Get problem dimension*/ + int extrapolationcase=GetExtrapolationCase(element); + switch(extrapolationcase){ + case 0: + if(!element->IsOnBase()) return NULL; + workelement = element->SpawnBasalElement(); + break; + case 1: case 2: case 3: workelement=element; break; + } + /*Fetch number of nodes */ - int numnodes = basalelement->GetNumberOfNodes(); + int numnodes = workelement->GetNumberOfNodes(); /*Initialize Element vector*/ - ElementVector* pe = basalelement->NewElementVector(); - for(i=0;iNewElementVector(); + for(int i=0;ivalues[i]=0.; - basalelement->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;}; return pe; }/*}}}*/ -void ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -256,13 +261,12 @@ /*Build B: */ - for(int i=0;i(basis); }/*}}}*/ -void ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system @@ -279,12 +283,11 @@ /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ - for(int i=0;iFindParam(&extrapolationvariable, ExtrapolationVariableEnum); + switch(extrapolationcase){ + case 0: + element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable); + break; + case 1: + element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable); + break; + case 2: + element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable); + break; + case 3: + element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable); + break; + } +}/*}}}*/ +int ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/ + + /* Get case of extrapolation, depending on domain quality, and extrapolation variable */ int domaintype, extrapolationvariable; + int extrapolationcase; element->FindParam(&domaintype,DomainTypeEnum); - element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum); switch(domaintype){ - case Domain2DhorizontalEnum: - element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable); - break; - case Domain3DEnum: - element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable); - break; - default: _error_("mesh "<FindParam(&extrapolationvariable, ExtrapolationVariableEnum); + if(extrapolationvariable==ThicknessEnum) extrapolationcase=2; // scalar fields that are constant along z-axis + else extrapolationcase=3; // scalar fields that vary along z-axis + break; + } + return extrapolationcase; }/*}}}*/ void ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 20116) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 20117) @@ -26,9 +26,10 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim); + void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); + int GetExtrapolationCase(Element* element); void SetConstraintsOnIce(Element* element); void UpdateConstraints(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 20116) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 20117) @@ -95,6 +95,5 @@ analysis_enum==HydrologyDCInefficientAnalysisEnum || analysis_enum==HydrologyDCEfficientAnalysisEnum || - analysis_enum==LevelsetAnalysisEnum || - analysis_enum==ExtrapolationAnalysisEnum + analysis_enum==LevelsetAnalysisEnum ){ if(iomodel->domaintype!=Domain2DhorizontalEnum){ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 20116) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 20117) @@ -24,4 +24,5 @@ void surfaceslope_core(FemModel* femmodel); void levelsetfunctionslope_core(FemModel* femmodel); +void levelset_core(FemModel* femmodel); void bedslope_core(FemModel* femmodel); void meshdeformation_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/levelset_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/levelset_core.cpp (revision 20117) +++ /issm/trunk-jpl/src/c/cores/levelset_core.cpp (revision 20117) @@ -0,0 +1,78 @@ +/*!\file: levelset_core.cpp + * \brief: levelset-module to update the ice domain + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void levelset_core(FemModel* femmodel){ + + /* intermediaries */ + bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,iscalving; + int domaintype, num_extrapol_vars, index; + int* extrapol_vars=NULL; + Analysis *analysis=NULL; + + /* recover parameters */ + femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum); + femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum); + femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); + femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum); + femmodel->parameters->FindParam(&iscalving,TransientIscalvingEnum); + if(isthermal && domaintype==Domain3DEnum) femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + + if(!islevelset) return; + + /* start the work from here */ + if(iscalving) Calvingx(femmodel); + if(VerboseSolution()) _printf0_(" computing levelset transport\n"); + /* smoothen slope of lsf for computation of normal on ice domain*/ + levelsetfunctionslope_core(femmodel); + + /* determine variables for extrapolation */ + num_extrapol_vars=0; + if(isstressbalance) num_extrapol_vars+=2; + if(ismasstransport) num_extrapol_vars+=1; + if(isthermal && domaintype==Domain3DEnum) num_extrapol_vars+=1; + extrapol_vars=xNew(num_extrapol_vars); + index=0; + if(isstressbalance){ + extrapol_vars[index]=VxEnum; index++; + extrapol_vars[index]=VyEnum; index++; + } + if(ismasstransport){extrapol_vars[index]=ThicknessEnum; index++;} + if(isthermal && domaintype==Domain3DEnum){ + if(isenthalpy){extrapol_vars[index]=EnthalpyEnum;} + else{extrapol_vars[index]=TemperatureEnum;} + index++; + } + + /* extrapolate */ + analysis = new ExtrapolationAnalysis(); + for(int iv=0;ivparameters->SetParam(extrapol_vars[iv],ExtrapolationVariableEnum); + analysis->Core(femmodel); + } + xDelete(extrapol_vars); + delete analysis; + + /* solve level set equation */ + analysis = new LevelsetAnalysis(); + analysis->Core(femmodel); + delete analysis; + + /* update vertices included for next calculation */ + GetMaskOfIceVerticesLSMx(femmodel); + + /* add computation domain mask to outputs */ + if(save_results){ + int outputs[1] = {IceMaskNodeActivationEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + } +} Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 20116) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 20117) @@ -102,33 +102,7 @@ if(isdamageevolution) damage_core(femmodel); - if(islevelset){ - if(iscalving) Calvingx(femmodel); - if(VerboseSolution()) _printf0_(" computing levelset transport\n"); - /* smoothen slope of lsf for computation of normal on ice domain*/ - levelsetfunctionslope_core(femmodel); + if(islevelset) levelset_core(femmodel); - /* extrapolate velocities onto domain with no ice */ - Analysis* extanalysis = new ExtrapolationAnalysis(); - const int nvars=3; - int vars[nvars] = {VxEnum, VyEnum, ThicknessEnum}; - for(int iv=0;ivparameters->SetParam(vars[iv],ExtrapolationVariableEnum); - extanalysis->Core(femmodel); - } - delete extanalysis; - - /* solve level set equation */ - analysis = new LevelsetAnalysis(); - analysis->Core(femmodel); - delete analysis; - - /* update vertices included for next calculation */ - GetMaskOfIceVerticesLSMx(femmodel); - - /* add computation domain mask to outputs */ - int outputs[1] = {IceMaskNodeActivationEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); - } - + /* from here on, prepare geometry for next time step*/ if(issmb)smb_core(femmodel); @@ -155,5 +129,5 @@ } - /*Calculate new Basal melting on Floating ice*/ + /*Calculate new basal melting on floating ice*/ FloatingiceMeltingRatex(femmodel);