Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17097) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17098) @@ -7,25 +7,52 @@ int LevelsetAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ - _error_("not implemented yet"); + return 1; } /*}}}*/ void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ - _error_("not implemented yet"); + // parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum)); } /*}}}*/ void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ - _error_("not implemented yet"); + int stabilization,finiteelement; + + /*Finite element type*/ + finiteelement = P1Enum; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + } /*}}}*/ void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - _error_("not implemented yet"); + int finiteelement=P1Enum; + ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement); } /*}}}*/ void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediary*/ + + /*Fetch parameters: */ + } /*}}}*/ void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediary*/ + + /*Fetch parameters: */ + }/*}}}*/ @@ -35,14 +62,132 @@ }/*}}}*/ ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/ - _error_("not implemented yet"); + /*Default, return NULL*/ + return NULL; }/*}}}*/ ElementMatrix* LevelsetAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + /* Jacobian required for the Newton solver */ _error_("not implemented yet"); }/*}}}*/ ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + /*Intermediaries */ + + int dim = 2; // solve for LSF in horizontal plane only + IssmDouble kappa; + IssmDouble Jdet, dt, D_scalar; + IssmDouble u,v,um,vm,ub,vb,vx,vy; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(dim*numnodes); + IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble D[dim][dim]; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input); + Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + D_scalar=gauss->weight*Jdet; + + /* Transient */ + if(dt!=0.){ + element->NodalFunctions(basis,gauss); + TripleMultiply(basis,numnodes,1,0, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + D_scalar*=dt; + } + + /* Advection */ + GetB(B,element,xyz_list,gauss); + GetBprime(Bprime,element,xyz_list,gauss); + vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; + vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; + ub=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here + vb=0.; + D[0][0]=D_scalar*(vx+ub); + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=D_scalar*(vy+vb); + TripleMultiply(B,dim,numnodes,1, + &D[0][0],dim,dim,0, + Bprime,dim,numnodes,0, + &Ke->values[0],1); + + /* Artificial Diffusion */ + kappa=0.; //FIXME: insert suitable value for kappa + GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed? + D[0][0]=D_scalar*kappa; + D[1][1]=D_scalar*kappa; + TripleMultiply(Bprime,dim,numnodes,1, + &D[0][0],dim,dim,0, + Bprime,dim,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; }/*}}}*/ ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediaries */ + int i, ig; + IssmDouble Jdet,dt; + IssmDouble phi; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector*/ + ElementVector* pe = element->NewElementVector(); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + if(dt!=0.){ + /*Initialize basis vector*/ + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* levelset_input = element->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + levelset_input->GetInputValue(&phi,gauss); + for(i=0;ivalues[i]+=Jdet*gauss->weight*phi*basis[i]; + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + } + else{for(i=0;ivalues[i]=0.;} + return pe; }/*}}}*/ void LevelsetAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ @@ -52,2 +197,58 @@ _error_("not implemented yet"); }/*}}}*/ +void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*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 + * by: + * Bi=[ N ] + * [ N ] + * where N is the finiteelement function for node i. + * + * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*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 + * by: + * Bi_prime=[ dN/dx ] + * [ dN/dy ] + * where N is the finiteelement function for node i. + * + * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(2*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B': */ + for(int i=0;i(dbasis); + +}/*}}}*/ + Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17097) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17098) @@ -28,4 +28,7 @@ void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + }; #endif Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17097) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17098) @@ -273,5 +273,4 @@ TransientIsgiaEnum, TransientIsdamageEnum, - TransientIslevelsetEnum, TransientNumRequestedOutputsEnum, TransientRequestedOutputsEnum, @@ -321,5 +320,4 @@ DamageEvolutionSolutionEnum, DamageEvolutionAnalysisEnum, - LevelsetAnalysisEnum, StressbalanceAnalysisEnum, StressbalanceSIAAnalysisEnum, @@ -674,4 +672,8 @@ LliboutryDuvalEnum, /*}}}*/ + /*Levelset related enums (will be moved to appropriate place when finished){{{*/ + LevelsetAnalysisEnum, + TransientIslevelsetEnum, + /*}}}*/ MaximumNumberOfDefinitionsEnum }; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17097) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17098) @@ -279,5 +279,4 @@ case TransientIsgiaEnum : return "TransientIsgia"; case TransientIsdamageEnum : return "TransientIsdamage"; - case TransientIslevelsetEnum : return "TransientIslevelset"; case TransientNumRequestedOutputsEnum : return "TransientNumRequestedOutputs"; case TransientRequestedOutputsEnum : return "TransientRequestedOutputs"; @@ -322,5 +321,4 @@ case DamageEvolutionSolutionEnum : return "DamageEvolutionSolution"; case DamageEvolutionAnalysisEnum : return "DamageEvolutionAnalysis"; - case LevelsetAnalysisEnum : return "LevelsetAnalysis"; case StressbalanceAnalysisEnum : return "StressbalanceAnalysis"; case StressbalanceSIAAnalysisEnum : return "StressbalanceSIAAnalysis"; @@ -634,4 +632,6 @@ case ArrheniusEnum : return "Arrhenius"; case LliboutryDuvalEnum : return "LliboutryDuval"; + case LevelsetAnalysisEnum : return "LevelsetAnalysis"; + case TransientIslevelsetEnum : return "TransientIslevelset"; case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions"; default : return "unknown"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17097) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17098) @@ -285,5 +285,4 @@ else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; else if (strcmp(name,"TransientIsdamage")==0) return TransientIsdamageEnum; - else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum; else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; @@ -328,5 +327,4 @@ else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; - else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; @@ -383,10 +381,10 @@ else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; + else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + else if (strcmp(name,"Contour")==0) return ContourEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BoolParam")==0) return BoolParamEnum; - else if (strcmp(name,"Contour")==0) return ContourEnum; - else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; + if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; @@ -506,10 +504,10 @@ else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; + else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; - else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; - else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; + if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; @@ -629,10 +627,10 @@ else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; + else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; + else if (strcmp(name,"XY")==0) return XYEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; - else if (strcmp(name,"XY")==0) return XYEnum; - else if (strcmp(name,"XYZ")==0) return XYZEnum; + if (strcmp(name,"XYZ")==0) return XYZEnum; else if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; @@ -649,4 +647,6 @@ else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; + else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum; else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum; else stage=7; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17097) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17098) @@ -271,5 +271,4 @@ def TransientIsgiaEnum(): return StringToEnum("TransientIsgia")[0] def TransientIsdamageEnum(): return StringToEnum("TransientIsdamage")[0] -def TransientIslevelsetEnum(): return StringToEnum("TransientIslevelset")[0] def TransientNumRequestedOutputsEnum(): return StringToEnum("TransientNumRequestedOutputs")[0] def TransientRequestedOutputsEnum(): return StringToEnum("TransientRequestedOutputs")[0] @@ -314,5 +313,4 @@ def DamageEvolutionSolutionEnum(): return StringToEnum("DamageEvolutionSolution")[0] def DamageEvolutionAnalysisEnum(): return StringToEnum("DamageEvolutionAnalysis")[0] -def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0] def StressbalanceAnalysisEnum(): return StringToEnum("StressbalanceAnalysis")[0] def StressbalanceSIAAnalysisEnum(): return StringToEnum("StressbalanceSIAAnalysis")[0] @@ -626,3 +624,5 @@ def ArrheniusEnum(): return StringToEnum("Arrhenius")[0] def LliboutryDuvalEnum(): return StringToEnum("LliboutryDuval")[0] +def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0] +def TransientIslevelsetEnum(): return StringToEnum("TransientIslevelset")[0] def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]