Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17138) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17139) @@ -370,5 +370,7 @@ ./cores/transient_core.cpp\ ./analyses/LevelsetAnalysis.h\ - ./analyses/LevelsetAnalysis.cpp + ./analyses/LevelsetAnalysis.cpp\ + ./analyses/ExtrapolationAnalysis.h\ + ./analyses/ExtrapolationAnalysis.cpp #}}} #Steadystate sources {{{ Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 17139) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 17139) @@ -0,0 +1,232 @@ +#include "./ExtrapolationAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +int ExtrapolationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ + return 1; +} +/*}}}*/ +void ExtrapolationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + _error_("not implemented yet"); +} +/*}}}*/ +void ExtrapolationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + 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,ExtrapolationVariableEnum);//FIXME: is this the correct way? +} +/*}}}*/ +void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + int finiteelement=P1Enum; + ::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement); +} +/*}}}*/ +void ExtrapolationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + _error_("not implemented yet"); + +} +/*}}}*/ +void ExtrapolationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + _error_("not implemented yet"); + +}/*}}}*/ + +/*Finite element Analysis*/ +void ExtrapolationAnalysis::Core(FemModel* femmodel){/*{{{*/ + + /*activate formulation: */ + femmodel->SetCurrentConfiguration(ExtrapolationAnalysisEnum); + + if(VerboseSolution()) _printf0_("extrapolation: call computational core:\n"); + solutionsequence_linear(femmodel); + +}/*}}}*/ +ElementVector* ExtrapolationAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* ExtrapolationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + /* Jacobian required for the Newton solver */ + _error_("not implemented yet"); +}/*}}}*/ +ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Intermediaries */ + const int dim = 2; + int i,row,col,stabilization; + IssmDouble Jdet,D_scalar,h; + IssmDouble dlevelset[dim],normal[dim]; + IssmDouble norm_dlevelset; + 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* B = xNew(dim*numnodes); + IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble D[dim][dim]; + + /*Retrieve all inputs and parameters*/ + Input* levelset_input=element->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input); + element->GetVerticesCoordinates(&xyz_list); + h = element->CharacteristicLength(); + + /* 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); + GetB(B,element,xyz_list,gauss); + GetBprime(Bprime,element,xyz_list,gauss); + + /* Get normal on node */ + levelset_input->GetInputDerivativeValue(&dlevelset[0],xyz_list,gauss); + norm_dlevelset=0.; + for(i=0;iweight*Jdet; + + for(row=0;rowvalues[0],1); + + stabilization=1; + if (stabilization==0){/* no stabilization, do nothing*/} + else if(stabilization==1){ + /*Streamline upwinding*/ + for(row=0;rowvalues[0],1); + } + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; + + }/*}}}*/ +ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/ + + /*Intermediaries */ + int i; + + /*Fetch number of nodes */ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector*/ + ElementVector* pe = element->NewElementVector(); + for(i=0;ivalues[i]=0.; + return pe; +}/*}}}*/ +void ExtrapolationAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + int meshtype, extrapolationvariable; + element->FindParam(&meshtype,MeshTypeEnum); + element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable); + break; + case Mesh3DEnum: + element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void ExtrapolationAnalysis::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/ExtrapolationAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 17139) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 17139) @@ -0,0 +1,34 @@ +/*! \file ExtrapolationAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _ExtrapolationAnalysis_ +#define _ExtrapolationAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class ExtrapolationAnalysis: public Analysis{ + + public: + /*Model processing*/ + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + 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/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 17138) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 17139) @@ -17,4 +17,5 @@ #include "./ExtrudeFromBaseAnalysis.h" #include "./ExtrudeFromTopAnalysis.h" +#include "./ExtrapolationAnalysis.h" #include "./FreeSurfaceBaseAnalysis.h" #include "./FreeSurfaceTopAnalysis.h" Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17138) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17139) @@ -114,5 +114,5 @@ if(isthermal && meshtype==Mesh3DEnum){ - if(VerboseSolution()) _printf0_(" computing temperatures\n"); + if(VerboseSolution()) _printf0_(" computing thermal regime\n"); #ifdef _HAVE_THERMAL_ thermal_core(femmodel); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17138) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17139) @@ -676,4 +676,6 @@ LevelsetAnalysisEnum, TransientIslevelsetEnum, + ExtrapolationAnalysisEnum, + ExtrapolationVariableEnum, /*}}}*/ MaximumNumberOfDefinitionsEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17138) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17139) @@ -635,4 +635,6 @@ case LevelsetAnalysisEnum : return "LevelsetAnalysis"; case TransientIslevelsetEnum : return "TransientIslevelset"; + case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis"; + case ExtrapolationVariableEnum : return "ExtrapolationVariable"; 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 17138) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17139) @@ -650,4 +650,6 @@ else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum; + else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; + else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; 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 17138) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17139) @@ -627,3 +627,5 @@ def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0] def TransientIslevelsetEnum(): return StringToEnum("TransientIslevelset")[0] +def ExtrapolationAnalysisEnum(): return StringToEnum("ExtrapolationAnalysis")[0] +def ExtrapolationVariableEnum(): return StringToEnum("ExtrapolationVariable")[0] def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]