Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17471) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17472) @@ -446,4 +446,6 @@ ./classes/Inputs/PentaInput.h\ ./classes/Inputs/PentaInput.cpp\ + ./classes/Inputs/TetraInput.h\ + ./classes/Inputs/TetraInput.cpp\ #}}} #DAKOTA sources {{{ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17472) @@ -20,4 +20,5 @@ switch(meshtype){ case Mesh3DEnum: numdofs=2; break; + case Mesh3DtetrasEnum: numdofs=2; break; case Mesh2DhorizontalEnum: numdofs=2; break; case Mesh2DverticalEnum: numdofs=1; break; @@ -29,4 +30,5 @@ switch(meshtype){ case Mesh3DEnum: numdofs=2; break; + case Mesh3DtetrasEnum: numdofs=2; break; case Mesh2DverticalEnum: numdofs=1; break; default: _error_("mesh type not supported yet"); @@ -37,4 +39,5 @@ switch(meshtype){ case Mesh3DEnum: numdofs=3; break; + case Mesh3DtetrasEnum: numdofs=3; break; case Mesh2DverticalEnum: numdofs=2; break; default: _error_("mesh type not supported yet"); @@ -45,4 +48,5 @@ switch(meshtype){ case Mesh3DEnum: numdofs=4; break; + case Mesh3DtetrasEnum: numdofs=4; break; case Mesh2DverticalEnum: numdofs=3; break; default: _error_("mesh type not supported yet"); @@ -215,6 +219,14 @@ if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); } + if(iomodel->meshtype==Mesh3DtetrasEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,LoadingforceZEnum); + iomodel->FetchDataToInput(elements,VzEnum,0.); + if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); + } if(iomodel->meshtype==Mesh2DverticalEnum){ - iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); } if(isFS){ @@ -963,4 +975,5 @@ case Mesh2DverticalEnum: dim = 2; dofpernode = 1; break; case Mesh3DEnum: dim = 3; dofpernode = 2; break; + case Mesh3DtetrasEnum: dim = 3; dofpernode = 2; break; default: _error_("mesh "<FindParam(&meshtype,MeshTypeEnum); switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; case Mesh3DEnum: dim = 3; break; - case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DtetrasEnum: dim = 3; break; default: _error_("mesh "<GetNumberOfNodes(); + for(int i=0;iSid(); + } +} +/*}}}*/ +void Element::GetNodesLidList(int* lidlist){/*{{{*/ + + _assert_(lidlist); + _assert_(nodes); + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iLid(); + } +} +/*}}}*/ +void Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/ + + int numvertices = this->GetNumberOfVertices(); + IssmDouble* xyz_list = xNew(numvertices*3); + ::GetVerticesCoordinates(xyz_list,this->vertices,numvertices); + + *pxyz_list = xyz_list; + +}/*}}}*/ void Element::GetVerticesSidList(int* sidlist){/*{{{*/ @@ -439,4 +468,103 @@ int numvertices = this->GetNumberOfVertices(); for(int i=0;ivertices[i]->Connectivity(); +} +/*}}}*/ +void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/ + this->inputs->ChangeEnum(original_enum,new_enum); +} +/*}}}*/ +void Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/ + + /*Intermediaries*/ + int i,t,row; + IssmDouble time; + + /*Branch on type of vector: nodal or elementary: */ + if(vector_type==1){ //nodal vector + + int numvertices = this->GetNumberOfVertices(); + int *vertexids = xNew(numvertices); + IssmDouble *values = xNew(numvertices); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(i=0;i(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + } + + /*Are we in transient or static? */ + if(M==iomodel->numberofvertices){ + for(i=0;iAddInput(vector_enum,values,P1Enum); + } + else if(M==iomodel->numberofvertices+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + + xDelete(values); + xDelete(vertexids); + } + else if(vector_type==2){ //element vector + /*Are we in transient or static? */ + if(M==iomodel->numberofelements){ + if (code==5){ //boolean + this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==6){ //integer + this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==7){ //IssmDouble + this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()])); + } + else _error_("could not recognize nature of vector from code " << code); + } + else { + _error_("transient element inputs not supported yet!"); + } + } + else{ + _error_("Cannot add input for vector type " << vector_type << " (not supported)"); + } +}/*}}}*/ +void Element::InputUpdateFromConstant(int constant, int name){/*{{{*/ + + /*Check that name is an element input*/ + if(!IsInput(name)) return; + + /*update input*/ + this->inputs->AddInput(new IntInput(name,constant)); +} +/*}}}*/ +void Element::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/ + + /*Check that name is an element input*/ + if(!IsInput(name)) return; + + /*update input*/ + this->inputs->AddInput(new DoubleInput(name,constant)); +} +/*}}}*/ +void Element::InputUpdateFromConstant(bool constant, int name){/*{{{*/ + + /*Check that name is an element input*/ + if(!IsInput(name)) return; + + /*update input*/ + this->inputs->AddInput(new BoolInput(name,constant)); } /*}}}*/ @@ -535,4 +663,16 @@ } /*}}}*/ +ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/ + return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/ + return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/ + return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum); +} +/*}}}*/ void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17472) @@ -75,9 +75,20 @@ void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type); IssmDouble GetMaterialParameter(int enum_in); + void GetNodesSidList(int* sidlist); + void GetNodesLidList(int* lidlist); void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); + void GetVerticesCoordinates(IssmDouble** xyz_list); void GetVerticesSidList(int* sidlist); void GetVerticesConnectivityList(int* connectivitylist); + void InputChangeName(int enum_type,int enum_type_old); + void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); + void InputUpdateFromConstant(IssmDouble constant, int name); + void InputUpdateFromConstant(int constant, int name); + void InputUpdateFromConstant(bool constant, int name); bool IsInput(int name); bool IsFloating(); + ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); + ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); + ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum); void ResultToVector(Vector* vector,int output_enum); @@ -161,6 +172,4 @@ virtual int GetNumberOfNodesPressure(void)=0; virtual int GetNumberOfVertices(void)=0; - virtual void GetNodesSidList(int* sidlist)=0; - virtual void GetNodesLidList(int* sidlist)=0; virtual int Id()=0; @@ -174,5 +183,4 @@ virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; virtual Node* GetNode(int node_number)=0; - virtual void GetVerticesCoordinates(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0; @@ -184,5 +192,4 @@ virtual IssmDouble SurfaceArea(void)=0; virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type)=0; - virtual void InputChangeName(int enum_type,int enum_type_old)=0; virtual void ComputeBasalStress(Vector* sigma_b)=0; virtual void ComputeStrainRate(Vector* eps)=0; @@ -191,5 +198,4 @@ virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0; virtual void InputDuplicate(int original_enum,int new_enum)=0; - virtual void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0; virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0; virtual void InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum)=0; @@ -206,7 +212,5 @@ virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0; virtual Gauss* NewGaussTop(int order)=0; - virtual ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum)=0; - virtual ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum)=0; - virtual ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum)=0; + virtual void InputScale(int enum_type,IssmDouble scale_factor)=0; virtual void GetVectorFromInputs(Vector* vector, int name_enum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17472) @@ -28,5 +28,5 @@ Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels) :PentaRef(nummodels) - ,ElementHook(nummodels,index+1,6,iomodel){ + ,ElementHook(nummodels,index+1,NUMVERTICES,iomodel){ int penta_elements_ids[2]; @@ -843,28 +843,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetNodesSidList{{{*/ -void Penta::GetNodesSidList(int* sidlist){ - - _assert_(sidlist); - _assert_(nodes); - int numnodes = this->NumberofNodes(); - - for(int i=0;iSid(); - } -} -/*}}}*/ -/*FUNCTION Penta::GetNodesLidList{{{*/ -void Penta::GetNodesLidList(int* lidlist){ - - _assert_(lidlist); - _assert_(nodes); - int numnodes = this->NumberofNodes(); - - for(int i=0;iLid(); - } -} -/*}}}*/ /*FUNCTION Penta::GetNumberOfNodes;{{{*/ int Penta::GetNumberOfNodes(void){ @@ -907,14 +883,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/ -void Penta::GetVerticesCoordinates(IssmDouble** pxyz_list){ - - IssmDouble* xyz_list = xNew(NUMVERTICES*3); - ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES); - - /*Assign output pointer*/ - *pxyz_list = xyz_list; - -}/*}}}*/ /*FUNCTION Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/ void Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){ @@ -953,5 +919,5 @@ cross(normal,AB,AC); - norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0)); + norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); for(int i=0;i<3;i++) normal[i]=normal[i]/norm; @@ -1206,81 +1172,4 @@ int Penta::Id(void){ return id; -} -/*}}}*/ -/*FUNCTION Penta::InputChangeName{{{*/ -void Penta::InputChangeName(int original_enum,int new_enum){ - - /*Call inputs method*/ - this->inputs->ChangeEnum(original_enum,new_enum); -} -/*}}}*/ -/*FUNCTION Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ -void Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){ - - /*Intermediaries*/ - int i,t; - int penta_vertex_ids[6]; - int row; - IssmDouble nodeinputs[6]; - IssmDouble time; - TransientInput *transientinput = NULL; - - /*Branch on type of vector: nodal or elementary: */ - if(vector_type==1){ //nodal vector - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<6;i++){ - _assert_(iomodel->elements); - penta_vertex_ids[i]=iomodel->elements[6*this->sid+i]; //ids for vertices are in the elements array from Matlab - } - - /*Are we in transient or static? */ - if(M==iomodel->numberofvertices){ - - /*create input values: */ - for(i=0;i<6;i++)nodeinputs[i]=(IssmDouble)vector[penta_vertex_ids[i]-1]; - - /*create static input: */ - this->inputs->AddInput(new PentaInput(vector_enum,nodeinputs,P1Enum)); - } - else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tAddTimeInput(new PentaInput(vector_enum,nodeinputs,P1Enum)); - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else _error_("nodal vector is either numberofvertices (" << iomodel->numberofvertices << "), or numberofvertices+1 long. Field provided is " << M << " long. Enum " << EnumToStringx(vector_enum)); - } - else if(vector_type==2){ //element vector - /*Are we in transient or static? */ - if(M==iomodel->numberofelements){ - - /*static mode: create an input out of the element value: */ - - if (code==5){ //boolean - this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->sid]))); - } - else if (code==6){ //integer - this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->sid]))); - } - else if (code==7){ //IssmDouble - this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->sid])); - } - else _error_("could not recognize nature of vector from code " << code); - } - else { - _error_("transient elementary inputs not supported yet!"); - } - } - else{ - _error_("Cannot add input for vector type " << vector_type << " (not supported)"); - } - } /*}}}*/ @@ -1452,41 +1341,13 @@ } /*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{*/ -void Penta::InputUpdateFromConstant(bool constant, int name){ - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new BoolInput(name,constant)); -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(IssmDouble value, int name);{{{*/ -void Penta::InputUpdateFromConstant(IssmDouble constant, int name){ - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new DoubleInput(name,constant)); -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{*/ -void Penta::InputUpdateFromConstant(int constant, int name){ - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new IntInput(name,constant)); -} -/*}}}*/ /*FUNCTION Penta::InputUpdateFromIoModel {{{*/ void Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ /*Intermediaries*/ - IssmInt i,j; - int penta_vertex_ids[6]; - IssmDouble nodeinputs[6]; - IssmDouble cmmininputs[6]; - IssmDouble cmmaxinputs[6]; + int i,j; + int penta_vertex_ids[NUMVERTICES]; + IssmDouble nodeinputs[NUMVERTICES]; + IssmDouble cmmininputs[NUMVERTICES]; + IssmDouble cmmaxinputs[NUMVERTICES]; IssmDouble yts; @@ -1500,12 +1361,8 @@ if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum); - /*Checks if debuging*/ - /*{{{*/ + /*Recover vertices ids needed to initialize inputs*/ _assert_(iomodel->elements); - /*}}}*/ - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<6;i++){ - penta_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + for(i=0;ielements[NUMVERTICES*index+i]; //ids for vertices are in the elements array from Matlab } @@ -1516,7 +1373,7 @@ case BalancethicknessThickeningRateEnum: if (iomodel->Data(BalancethicknessThickeningRateEnum)){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[penta_vertex_ids[j]-1]/yts; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(BalancethicknessThickeningRateEnum)[penta_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1524,7 +1381,7 @@ case VxEnum: if (iomodel->Data(VxEnum)){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(VxEnum)[penta_vertex_ids[j]-1]/yts; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(VxEnum)[penta_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1532,7 +1389,7 @@ case VyEnum: if (iomodel->Data(VyEnum)){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(VyEnum)[penta_vertex_ids[j]-1]/yts; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(VyEnum)[penta_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1540,7 +1397,7 @@ case FrictionCoefficientEnum: if (iomodel->Data(FrictionCoefficientEnum)){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[penta_vertex_ids[j]-1]; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(FrictionCoefficientEnum)[penta_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1548,7 +1405,7 @@ case MaterialsRheologyBbarEnum: if(iomodel->Data(MaterialsRheologyBEnum)){ - for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1]; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1556,7 +1413,7 @@ case DamageDbarEnum: if(iomodel->Data(DamageDEnum)){ - for(j=0;j<6;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[penta_vertex_ids[j]-1]; - for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; - for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(DamageDEnum)[penta_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } @@ -1579,5 +1436,5 @@ DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum); for(i=0;iData(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i]; + for(j=0;jData(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i]; datasetinput->AddInput(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast(iomodel->Data(InversionCostFunctionsEnum)[i])); } @@ -1887,19 +1744,4 @@ Gauss* Penta::NewGaussTop(int order){ return new GaussPenta(3,4,5,order); -} -/*}}}*/ -/*FUNCTION Penta::NewElementVector{{{*/ -ElementVector* Penta::NewElementVector(int approximation_enum){ - return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -/*FUNCTION Penta::NewElementMatrix{{{*/ -ElementMatrix* Penta::NewElementMatrix(int approximation_enum){ - return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -/*FUNCTION Penta::NewElementMatrixCoupling{{{*/ -ElementMatrix* Penta::NewElementMatrixCoupling(int number_nodes,int approximation_enum){ - return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum); } /*}}}*/ @@ -2212,7 +2054,6 @@ void Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){ + /*go into parameters and get the analysis_counter: */ int analysis_counter; - - /*go into parameters and get the analysis_counter: */ parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); @@ -2220,7 +2061,8 @@ this->element_type=this->element_type_list[analysis_counter]; - /*Pick up nodes */ - if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + /*Pick up nodes*/ + if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17472) @@ -50,7 +50,4 @@ /*}}}*/ /*Update virtual functions definitions: {{{*/ - void InputUpdateFromConstant(bool constant, int name); - void InputUpdateFromConstant(IssmDouble constant, int name); - void InputUpdateFromConstant(int constant, int name); void InputUpdateFromVector(IssmDouble* vector, int name, int type); #ifdef _HAVE_DAKOTA_ @@ -85,6 +82,4 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); int GetNodeIndex(Node* node); - void GetNodesSidList(int* sidlist); - void GetNodesLidList(int* lidlist); int GetNumberOfNodes(void); int GetNumberOfNodesPressure(void); @@ -96,10 +91,8 @@ IssmDouble GetZcoord(Gauss* gauss); void GetVectorFromInputs(Vector* vector,int name_enum); - void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); int Sid(); - void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputDuplicate(int original_enum,int new_enum); @@ -204,5 +197,4 @@ void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); Node* GetNode(int node_number); - void InputChangeName(int input_enum, int enum_type_old); void InputExtrude(int enum_type); void InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type); @@ -225,7 +217,4 @@ Gauss* NewGaussLine(int vertex1,int vertex2,int order); Gauss* NewGaussTop(int order); - ElementVector* NewElementVector(int approximation_enum); - ElementMatrix* NewElementMatrix(int approximation_enum); - ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17472) @@ -21,5 +21,5 @@ /*FUNCTION Seg::Seg(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/ Seg::Seg(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels) - :SegRef(nummodels),ElementHook(nummodels,index+1,2,iomodel){ + :SegRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){ /*id: */ @@ -263,14 +263,4 @@ } /*}}}*/ -/*FUNCTION Seg::NewElementVector{{{*/ -ElementVector* Seg::NewElementVector(int approximation_enum){ - return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -/*FUNCTION Seg::NewElementMatrix{{{*/ -ElementMatrix* Seg::NewElementMatrix(int approximation_enum){ - return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ /*FUNCTION Seg::NodalFunctions{{{*/ void Seg::NodalFunctions(IssmDouble* basis, Gauss* gauss){ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17472) @@ -53,7 +53,4 @@ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");}; #endif - void InputUpdateFromConstant(IssmDouble constant, int name){_error_("not implemented yet");}; - void InputUpdateFromConstant(int constant, int name){_error_("not implemented yet");}; - void InputUpdateFromConstant(bool constant, int name){_error_("not implemented yet");}; void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; /*}}}*/ @@ -80,6 +77,4 @@ Element* GetBasalElement(void){_error_("not implemented yet");}; int GetNodeIndex(Node* node){_error_("not implemented yet");}; - void GetNodesSidList(int* sidlist){_error_("not implemented yet");}; - void GetNodesLidList(int* lidlist){_error_("not implemented yet");}; int GetNumberOfNodes(void); int GetNumberOfNodesVelocity(void){_error_("not implemented yet");}; @@ -90,5 +85,4 @@ void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; int Sid(){_error_("not implemented yet");}; - void InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");}; bool IsOnBed(){_error_("not implemented yet");}; bool IsOnSurface(){_error_("not implemented yet");}; @@ -138,7 +132,4 @@ Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");}; - ElementVector* NewElementVector(int approximation_enum); - ElementMatrix* NewElementMatrix(int approximation_enum); - ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");}; int VertexConnectivity(int vertexindex){_error_("not implemented yet");}; void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; @@ -152,5 +143,4 @@ void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; void GetVectorFromInputs(Vector* vector, int name_enum){_error_("not implemented yet");}; - void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){_error_("not implemented yet");}; void InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");}; void InputDuplicate(int original_enum,int new_enum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17472) @@ -22,5 +22,5 @@ /*FUNCTION Tetra::Tetra(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/ Tetra::Tetra(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels) - :TetraRef(nummodels),ElementHook(nummodels,index+1,2,iomodel){ + :TetraRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){ /*id: */ @@ -127,4 +127,107 @@ /*}}}*/ +/*FUNCTION Tetra::AddInput{{{*/ +void Tetra::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){ + + /*Call inputs method*/ + _assert_(this->inputs); + this->inputs->AddInput(new TetraInput(input_enum,values,interpolation_enum)); +} +/*}}}*/ +/*FUNCTION Tetra::Configure {{{*/ +void Tetra::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){ + + int analysis_counter; + + /*go into parameters and get the analysis_counter: */ + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective + * datasets, using internal ids and offsets hidden in hooks: */ + if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); + this->hvertices->configure(verticesin); + this->hmaterial->configure(materialsin); + this->hmatpar->configure(materialsin); + + /*Now, go pick up the objects inside the hooks: */ + if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + this->vertices = (Vertex**)this->hvertices->deliverp(); + this->material = (Material*)this->hmaterial->delivers(); + this->matpar = (Matpar*)this->hmatpar->delivers(); + + /*point parameters to real dataset: */ + this->parameters=parametersin; + + /*get inputs configured too: */ + this->inputs->Configure(parameters); +} +/*}}}*/ +/*FUNCTION Tetra::FaceOnBedIndices{{{*/ +void Tetra::FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3){ + + IssmDouble values[NUMVERTICES]; + int indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}}; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonbedEnum); + + for(int i=0;i<4;i++){ + if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1. && values[indices[i][2]] == 1.){ + *pindex1 = indices[i][0]; + *pindex2 = indices[i][1]; + *pindex3 = indices[i][2]; + return; + } + } + + _error_("Could not find 3 vertices on bed"); +} +/*}}}*/ +/*FUNCTION Tetra::FaceOnSurfaceIndices{{{*/ +void Tetra::FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3){ + + IssmDouble values[NUMVERTICES]; + int indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}}; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum); + + for(int i=0;i<4;i++){ + if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1. && values[indices[i][2]] == 1.){ + *pindex1 = indices[i][0]; + *pindex2 = indices[i][1]; + *pindex3 = indices[i][2]; + return; + } + } + + _error_("Could not find 3 vertices on bed"); +} +/*}}}*/ +/*FUNCTION Tetra::FaceOnFranceIndices{{{*/ +void Tetra::FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3){ + + IssmDouble values[NUMVERTICES]; + int indices[4][3] = {{0,1,2},{0,1,3},{1,2,3},{2,0,3}}; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MaskIceLevelsetEnum); + + for(int i=0;i<4;i++){ + if(values[indices[i][0]] == 0. && values[indices[i][1]] == 0. && values[indices[i][2]] == 0.){ + *pindex1 = indices[i][0]; + *pindex2 = indices[i][1]; + *pindex3 = indices[i][2]; + return; + } + } + + _error_("Could not find 3 vertices on bed"); +} +/*}}}*/ /*FUNCTION Tetra::GetNumberOfNodes;{{{*/ int Tetra::GetNumberOfNodes(void){ @@ -137,23 +240,237 @@ } /*}}}*/ -/*FUNCTION Tetra::GetVerticesCoordinates(IssmDouble** pxyz_list) THIS ONE{{{*/ -void Tetra::GetVerticesCoordinates(IssmDouble** pxyz_list){ - - IssmDouble* xyz_list = xNew(NUMVERTICES*3); - ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES); +/*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/ +void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){ + + int indices[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Allocate Output*/ + IssmDouble* xyz_list_edge = xNew(3*3); + this->FaceOnBedIndices(&indices[0],&indices[1],&indices[2]); + for(int i=0;i<3;i++) for(int j=0;j<3;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; /*Assign output pointer*/ - *pxyz_list = xyz_list; + *pxyz_list = xyz_list_edge; }/*}}}*/ -/*FUNCTION Tetra::IsIceInElement {{{*/ +/*FUNCTION Tetra::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/ +void Tetra::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){ + + int indices[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Allocate Output*/ + IssmDouble* xyz_list_edge = xNew(3*3); + this->FaceOnSurfaceIndices(&indices[0],&indices[1],&indices[2]); + for(int i=0;i<3;i++) for(int j=0;j<3;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; + + /*Assign output pointer*/ + *pxyz_list = xyz_list_edge; + +}/*}}}*/ +/*FUNCTION Tetra::GetZcoord {{{*/ +IssmDouble Tetra::GetZcoord(Gauss* gauss){ + + int i; + IssmDouble z; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble z_list[NUMVERTICES]; + + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + for(i=0;iConstant(&yts,ConstantsYtsEnum); + iomodel->Constant(&control_analysis,InversionIscontrolEnum); + if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum); + if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(i=0;ielements[NUMVERTICES*index+i]; //ids for vertices are in the elements array from Matlab + } + + /*Control Inputs*/ + if (control_analysis && iomodel->Data(InversionControlParametersEnum)){ + for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i])){ + case BalancethicknessThickeningRateEnum: + if (iomodel->Data(BalancethicknessThickeningRateEnum)){ + for(j=0;jData(BalancethicknessThickeningRateEnum)[tetra_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case VxEnum: + if (iomodel->Data(VxEnum)){ + for(j=0;jData(VxEnum)[tetra_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case VyEnum: + if (iomodel->Data(VyEnum)){ + for(j=0;jData(VyEnum)[tetra_vertex_ids[j]-1]/yts; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; + this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case FrictionCoefficientEnum: + if (iomodel->Data(FrictionCoefficientEnum)){ + for(j=0;jData(FrictionCoefficientEnum)[tetra_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case MaterialsRheologyBbarEnum: + if(iomodel->Data(MaterialsRheologyBEnum)){ + for(j=0;jData(MaterialsRheologyBEnum)[tetra_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case DamageDbarEnum: + if(iomodel->Data(DamageDEnum)){ + for(j=0;jData(DamageDEnum)[tetra_vertex_ids[j]-1]; + for(j=0;jData(InversionMinParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[(tetra_vertex_ids[j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + default: + _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); + } + } + } + + /*Need to know the type of approximation for this element*/ + if(iomodel->Data(FlowequationElementEquationEnum)){ + this->inputs->AddInput(new IntInput(ApproximationEnum,reCast(iomodel->Data(FlowequationElementEquationEnum)[index]))); + } + + /*DatasetInputs*/ + if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)) { + + /*Create inputs and add to DataSetInput*/ + DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum); + for(i=0;iData(InversionCostFunctionsCoefficientsEnum)[(tetra_vertex_ids[j]-1)*num_responses+i]; + datasetinput->AddInput(new TetraInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast(iomodel->Data(InversionCostFunctionsEnum)[i])); + } + + /*Add datasetinput to element inputs*/ + this->inputs->AddInput(datasetinput); + } +} +/*}}}*/ +/*FUNCTION Tetra::IsIceInElement THIS ONE{{{*/ bool Tetra::IsIceInElement(){ - _error_("not implemented yet"); + /*Get levelset*/ + IssmDouble ls[NUMVERTICES]; + GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); + + /*If the level set on one of the nodes is <0, ice is present in this element*/ + for(int i=0;iEnum()==GaussTetraEnum); + this->GetJacobianDeterminantFace(pJdet,xyz_list,(GaussTetra*)gauss); + +} +/*}}}*/ +/*FUNCTION Tetra::JacobianDeterminantBase{{{*/ +void Tetra::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTetraEnum); + this->GetJacobianDeterminantFace(pJdet,xyz_list_base,(GaussTetra*)gauss); + +} +/*}}}*/ +/*FUNCTION Tetra::JacobianDeterminantTop{{{*/ +void Tetra::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTetraEnum); + this->GetJacobianDeterminantFace(pJdet,xyz_list_base,(GaussTetra*)gauss); } @@ -182,12 +516,26 @@ } /*}}}*/ -/*FUNCTION Tetra::NewElementVector THIS ONE{{{*/ -ElementVector* Tetra::NewElementVector(int approximation_enum){ - return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -/*FUNCTION Tetra::NewElementMatrix THIS ONE{{{*/ -ElementMatrix* Tetra::NewElementMatrix(int approximation_enum){ - return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum); +/*FUNCTION Tetra::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/ +Gauss* Tetra::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){ + /*FIXME: this is messed up, should provide indices and not xyz_list!*/ + int indices[3]; + this->FaceOnFrontIndices(&indices[0],&indices[1],&indices[2]); + return new GaussTetra(indices[0],indices[1],indices[2],max(order_horiz,order_vert)); +} +/*}}}*/ +/*FUNCTION Tetra::NewGaussBase(int order){{{*/ +Gauss* Tetra::NewGaussBase(int order){ + + int indices[3]; + this->FaceOnBedIndices(&indices[0],&indices[1],&indices[2]); + return new GaussTetra(indices[0],indices[1],indices[2],order); +} +/*}}}*/ +/*FUNCTION Tetra::NewGaussTop(int order){{{*/ +Gauss* Tetra::NewGaussTop(int order){ + + int indices[3]; + this->FaceOnSurfaceIndices(&indices[0],&indices[1],&indices[2]); + return new GaussTetra(indices[0],indices[1],indices[2],order); } /*}}}*/ @@ -209,7 +557,271 @@ /*}}}*/ /*FUNCTION Tetra::NormalSection{{{*/ -void Tetra::NormalSection(IssmDouble* normal,IssmDouble* xyz_list_front){ - - _error_("Not implemented yet"); -} -/*}}}*/ +void Tetra::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){ + + /*Build unit outward pointing vector*/ + IssmDouble AB[3]; + IssmDouble AC[3]; + IssmDouble norm; + + AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; + AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; + AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2]; + AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0]; + AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1]; + AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2]; + + cross(normal,AB,AC); + norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); + + for(int i=0;i<3;i++) normal[i]=normal[i]/norm; +} +/*}}}*/ +/*FUNCTION Tetra::ReduceMatrices{{{*/ +void Tetra::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){ + + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + if(pe){ + if(analysis_type==StressbalanceAnalysisEnum){ + if(this->element_type==MINIcondensedEnum){ + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ + //Do nothing, condensation already done in PVectorCoupling + } + else{ + _error_("Not implemented"); + } + } + else if(this->element_type==P1bubblecondensedEnum){ + _error_("Not implemented"); + } + } + } + + if(Ke){ + if(analysis_type==StressbalanceAnalysisEnum){ + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ + //Do nothing condensatino already done for Stokes part + } + else{ + if(this->element_type==MINIcondensedEnum){ + _error_("Not implemented"); + } + else if(this->element_type==P1bubblecondensedEnum){ + _error_("Not implemented"); + } + } + } + } +} +/*}}}*/ +/*FUNCTION Tetra::SetCurrentConfiguration {{{*/ +void Tetra::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){ + + /*go into parameters and get the analysis_counter: */ + int analysis_counter; + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Pick up nodes*/ + if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + +} +/*}}}*/ +/*FUNCTION Tetra::SetwiseNodeConnectivity THIS ONE (take from Penta.cpp){{{*/ +void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){ + + /*Intermediaries*/ + const int numnodes = this->NumberofNodes(); + + /*Output */ + int d_nz = 0; + int o_nz = 0; + + /*Loop over all nodes*/ + for(int i=0;inodes[i]->Lid()]){ + + /*flag current node so that no other element processes it*/ + flags[this->nodes[i]->Lid()]=true; + + int counter=0; + while(flagsindices[counter]>=0) counter++; + flagsindices[counter]=this->nodes[i]->Lid(); + + /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ + switch(set2_enum){ + case FsetEnum: + if(nodes[i]->indexing.fsize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + case GsetEnum: + if(nodes[i]->indexing.gsize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + case SsetEnum: + if(nodes[i]->indexing.ssize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + default: _error_("not supported"); + } + } + } + + /*Assign output pointers: */ + *pd_nz=d_nz; + *po_nz=o_nz; +} +/*}}}*/ +/*FUNCTION Tetra::Sid THIS ONE{{{*/ +int Tetra::Sid(){ + + return sid; + +} +/*}}}*/ +/*FUNCTION Tetra::Update {{{*/ +void Tetra::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ + + /*Intermediaries*/ + int i; + int tetra_vertex_ids[6]; + IssmDouble nodeinputs[6]; + IssmDouble yts; + bool dakota_analysis; + bool isFS; + int numnodes; + int* tetra_node_ids = NULL; + + /*Fetch parameters: */ + iomodel->Constant(&yts,ConstantsYtsEnum); + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + + /*Checks if debuging*/ + _assert_(iomodel->elements); + + /*Recover element type*/ + this->SetElementType(finiteelement_type,analysis_counter); + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<6;i++) tetra_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + + /*Recover nodes ids needed to initialize the node hook.*/ + switch(finiteelement_type){ + case P1Enum: + numnodes = 4; + tetra_node_ids = xNew(numnodes); + tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0]; + tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1]; + tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2]; + tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3]; + break; + default: + _error_("Finite element "<SetHookNodes(tetra_node_ids,numnodes,analysis_counter); this->nodes=NULL; + xDelete(tetra_node_ids); + + /*Fill with IoModel*/ + this->InputUpdateFromIoModel(index,iomodel); +} +/*}}}*/ +/*FUNCTION Tetra::ZeroLevelsetCoordinates{{{*/ +void Tetra::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){ + /*Compute portion of the element that is grounded*/ + + IssmDouble levelset[NUMVERTICES]; + IssmDouble* xyz_zero = xNew(3*3); + + /*Recover parameters and values*/ + GetInputListOnVertices(&levelset[0],levelsetenum); + + if(levelset[0]==0. && levelset[1]==0. && levelset[2]==0.){ + xyz_zero[3*0+0]=xyz_list[0*3+0]; + xyz_zero[3*0+1]=xyz_list[0*3+1]; + xyz_zero[3*0+2]=xyz_list[0*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[1*3+0]; + xyz_zero[3*1+1]=xyz_list[1*3+1]; + xyz_zero[3*1+2]=xyz_list[1*3+2]; + + /*New point 3*/ + xyz_zero[3*2+0]=xyz_list[2*3+0]; + xyz_zero[3*2+1]=xyz_list[2*3+1]; + xyz_zero[3*2+2]=xyz_list[2*3+2]; + } + else if(levelset[0]==0. && levelset[1]==0. && levelset[3]==0.){ + xyz_zero[3*0+0]=xyz_list[0*3+0]; + xyz_zero[3*0+1]=xyz_list[0*3+1]; + xyz_zero[3*0+2]=xyz_list[0*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[1*3+0]; + xyz_zero[3*1+1]=xyz_list[1*3+1]; + xyz_zero[3*1+2]=xyz_list[1*3+2]; + + /*New point 3*/ + xyz_zero[3*2+0]=xyz_list[3*3+0]; + xyz_zero[3*2+1]=xyz_list[3*3+1]; + xyz_zero[3*2+2]=xyz_list[3*3+2]; + } + else if(levelset[1]==0. && levelset[2]==0. && levelset[3]==0.){ + xyz_zero[3*0+0]=xyz_list[1*3+0]; + xyz_zero[3*0+1]=xyz_list[1*3+1]; + xyz_zero[3*0+2]=xyz_list[1*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[2*3+0]; + xyz_zero[3*1+1]=xyz_list[2*3+1]; + xyz_zero[3*1+2]=xyz_list[2*3+2]; + + /*New point 3*/ + xyz_zero[3*2+0]=xyz_list[3*3+0]; + xyz_zero[3*2+1]=xyz_list[3*3+1]; + xyz_zero[3*2+2]=xyz_list[3*3+2]; + } + else if(levelset[2]==0. && levelset[0]==0. && levelset[3]==0.){ + xyz_zero[3*0+0]=xyz_list[2*3+0]; + xyz_zero[3*0+1]=xyz_list[2*3+1]; + xyz_zero[3*0+2]=xyz_list[2*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[0*3+0]; + xyz_zero[3*1+1]=xyz_list[0*3+1]; + xyz_zero[3*1+2]=xyz_list[0*3+2]; + + /*New point 3*/ + xyz_zero[3*2+0]=xyz_list[3*3+0]; + xyz_zero[3*2+1]=xyz_list[3*3+1]; + xyz_zero[3*2+2]=xyz_list[3*3+2]; + } + else _error_("Case not covered"); + + /*Assign output pointer*/ + *pxyz_zero= xyz_zero; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17472) @@ -53,19 +53,16 @@ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type){_error_("not implemented yet");}; #endif - void InputUpdateFromConstant(IssmDouble constant, int name){_error_("not implemented yet");}; - void InputUpdateFromConstant(int constant, int name){_error_("not implemented yet");}; - void InputUpdateFromConstant(bool constant, int name){_error_("not implemented yet");}; - void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; + void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Element virtual functions definitions: {{{*/ void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");}; - void AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");}; + void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; void ComputeBasalStress(Vector* sigma_b){_error_("not implemented yet");}; void ComputeStrainRate(Vector* eps){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; - void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");}; + void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); + void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void Delta18oParameterization(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; @@ -74,4 +71,7 @@ IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; + void FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3); + void FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3); + void FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3); int FiniteElement(void); Element* GetUpperElement(void){_error_("not implemented yet");}; @@ -80,23 +80,21 @@ Element* GetBasalElement(void){_error_("not implemented yet");}; int GetNodeIndex(Node* node){_error_("not implemented yet");}; - void GetNodesSidList(int* sidlist){_error_("not implemented yet");}; - void GetNodesLidList(int* lidlist){_error_("not implemented yet");}; int GetNumberOfNodes(void); int GetNumberOfNodesVelocity(void){_error_("not implemented yet");}; int GetNumberOfNodesPressure(void){_error_("not implemented yet");}; int GetNumberOfVertices(void); - void GetVerticesCoordinates(IssmDouble** pxyz_list); - void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; - void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; - int Sid(){_error_("not implemented yet");}; - void InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");}; - bool IsOnBed(){_error_("not implemented yet");}; - bool IsOnSurface(){_error_("not implemented yet");}; + void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); + void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); + bool HasFaceOnBed(); + bool HasFaceOnSurface(); + int Sid(); + bool IsOnBed(); + bool IsOnSurface(); bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); - void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; - void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); + void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss); @@ -128,17 +126,14 @@ IssmDouble GetXcoord(Gauss* gauss){_error_("Not implemented");}; IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");}; - IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented yet");}; + IssmDouble GetZcoord(Gauss* gauss); int GetElementType(void){_error_("not implemented yet");}; Gauss* NewGauss(void); Gauss* NewGauss(int order); Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; - Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");}; + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");}; - Gauss* NewGaussBase(int order){_error_("not implemented yet");}; + Gauss* NewGaussBase(int order); Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; - Gauss* NewGaussTop(int order){_error_("not implemented yet");}; - ElementVector* NewElementVector(int approximation_enum); - ElementMatrix* NewElementMatrix(int approximation_enum); - ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");}; + Gauss* NewGaussTop(int order); int VertexConnectivity(int vertexindex){_error_("not implemented yet");}; void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; @@ -146,5 +141,5 @@ bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; bool IsIcefront(void); - void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; + void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; void GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");}; @@ -152,5 +147,4 @@ void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; void GetVectorFromInputs(Vector* vector, int name_enum){_error_("not implemented yet");}; - void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){_error_("not implemented yet");}; void InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");}; void InputDuplicate(int original_enum,int new_enum){_error_("not implemented yet");}; @@ -160,9 +154,9 @@ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; void ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");}; - void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");}; + void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; void SmbGradients(){_error_("not implemented yet");}; IssmDouble SurfaceArea(void){_error_("not implemented yet");}; - void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");}; + void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); IssmDouble TimeAdapt(){_error_("not implemented yet");}; void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp (revision 17472) @@ -329,4 +329,24 @@ } /*}}}*/ +/*FUNCTION TetraRef::GetJacobianDeterminantFace{{{*/ +void TetraRef::GetJacobianDeterminantFace(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTetra* gauss){ + /*The Jacobian determinant is constant over the element, discard the gaussian points. + * J is assumed to have been allocated of size NDOF2xNDOF2.*/ + + IssmDouble x1=xyz_list[3*0+0]; + IssmDouble y1=xyz_list[3*0+1]; + IssmDouble z1=xyz_list[3*0+2]; + IssmDouble x2=xyz_list[3*1+0]; + IssmDouble y2=xyz_list[3*1+1]; + IssmDouble z2=xyz_list[3*1+2]; + IssmDouble x3=xyz_list[3*2+0]; + IssmDouble y3=xyz_list[3*2+1]; + IssmDouble z3=xyz_list[3*2+2]; + + /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */ + *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5); + if(*Jdet<0) _error_("negative jacobian determinant!"); +} +/*}}}*/ /*FUNCTION TetraRef::GetJacobianInvert {{{*/ void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){ Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h (revision 17472) @@ -25,4 +25,5 @@ void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss); void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTetra* gauss); + void GetJacobianDeterminantFace(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTetra* gauss); void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss); void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17472) @@ -26,5 +26,5 @@ /*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{*/ Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels) - :TriaRef(nummodels),ElementHook(nummodels,index+1,3,iomodel){ + :TriaRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){ /*id: */ @@ -687,14 +687,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/ -void Tria::GetVerticesCoordinates(IssmDouble** pxyz_list){ - - IssmDouble* xyz_list = xNew(NUMVERTICES*3); - ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES); - - /*Assign output pointer*/ - *pxyz_list = xyz_list; - -}/*}}}*/ /*FUNCTION Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/ void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){ @@ -927,28 +917,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetNodesSidList{{{*/ -void Tria::GetNodesSidList(int* sidlist){ - - _assert_(sidlist); - _assert_(nodes); - int numnodes = this->NumberofNodes(); - - for(int i=0;iSid(); - } -} -/*}}}*/ -/*FUNCTION Tria::GetNodesLidList{{{*/ -void Tria::GetNodesLidList(int* lidlist){ - - _assert_(lidlist); - _assert_(nodes); - int numnodes = this->NumberofNodes(); - - for(int i=0;iLid(); - } -} -/*}}}*/ /*FUNCTION Tria::GetNumberOfNodes;{{{*/ int Tria::GetNumberOfNodes(void){ @@ -1086,11 +1052,4 @@ } /*}}}*/ -/*FUNCTION Tria::InputChangeName{{{*/ -void Tria::InputChangeName(int original_enum,int new_enum){ - - /*Call inputs method*/ - this->inputs->ChangeEnum(original_enum,new_enum); -} -/*}}}*/ /*FUNCTION Tria::InputScale{{{*/ void Tria::InputScale(int enum_type,IssmDouble scale_factor){ @@ -1104,34 +1063,4 @@ /*Scale: */ input->Scale(scale_factor); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{*/ -void Tria::InputUpdateFromConstant(int constant, int name){ - - /*Check that name is an element input*/ - if(!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new IntInput(name,constant)); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(IssmDouble value, int name);{{{*/ -void Tria::InputUpdateFromConstant(IssmDouble constant, int name){ - - /*Check that name is an element input*/ - if(!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new DoubleInput(name,constant)); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{*/ -void Tria::InputUpdateFromConstant(bool constant, int name){ - - /*Check that name is an element input*/ - if(!IsInput(name)) return; - - /*update input*/ - this->inputs->AddInput(new BoolInput(name,constant)); } /*}}}*/ @@ -1340,74 +1269,4 @@ } /*}}}*/ -/*FUNCTION Tria::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/ -void Tria::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){ - - /*Intermediaries*/ - int i,t; - int tria_vertex_ids[3]; - int row; - IssmDouble nodeinputs[3]; - IssmDouble time; - TransientInput* transientinput=NULL; - - /*Branch on type of vector: nodal or elementary: */ - if(vector_type==1){ //nodal vector - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<3;i++){ - _assert_(iomodel->elements); - tria_vertex_ids[i]=reCast(iomodel->elements[3*this->sid+i]); //ids for vertices are in the elements array from Matlab - } - - /*Are we in transient or static? */ - if(M==iomodel->numberofvertices){ - - /*create input values: */ - for(i=0;i<3;i++)nodeinputs[i]=vector[tria_vertex_ids[i]-1]; - - /*create static input: */ - this->inputs->AddInput(new TriaInput(vector_enum,nodeinputs,P1Enum)); - } - else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tAddTimeInput(new TriaInput(vector_enum,nodeinputs,P1Enum)); - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - } - else if(vector_type==2){ //element vector - /*Are we in transient or static? */ - if(M==iomodel->numberofelements){ - - /*static mode: create an input out of the element value: */ - - if (code==5){ //boolean - this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->sid]))); - } - else if (code==6){ //integer - this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->sid]))); - } - else if (code==7){ //IssmDouble - this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->sid])); - } - else _error_("could not recognize nature of vector from code " << code); - } - else { - _error_("transient element inputs not supported yet!"); - } - } - else{ - _error_("Cannot add input for vector type " << vector_type << " (not supported)"); - } - -} -/*}}}*/ /*FUNCTION Tria::IsOnBed {{{*/ bool Tria::IsOnBed(){ @@ -1681,14 +1540,4 @@ this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); return new GaussTria(indices[0],indices[1],order); -} -/*}}}*/ -/*FUNCTION Tria::NewElementVector{{{*/ -ElementVector* Tria::NewElementVector(int approximation_enum){ - return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -/*FUNCTION Tria::NewElementMatrix{{{*/ -ElementMatrix* Tria::NewElementMatrix(int approximation_enum){ - return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum); } /*}}}*/ @@ -2024,20 +1873,4 @@ } /*}}}*/ -/*FUNCTION Tria::SetCurrentConfiguration {{{*/ -void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){ - - /*go into parameters and get the analysis_counter: */ - int analysis_counter; - parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Get Element type*/ - this->element_type=this->element_type_list[analysis_counter]; - - /*Pick up nodes*/ - if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); - else this->nodes=NULL; - -} -/*}}}*/ /*FUNCTION Tria::SpawnSeg {{{*/ Seg* Tria::SpawnSeg(int index1,int index2){ @@ -2104,4 +1937,20 @@ _error_("not implemented yet"); } +} +/*}}}*/ +/*FUNCTION Tria::SetCurrentConfiguration {{{*/ +void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){ + + /*go into parameters and get the analysis_counter: */ + int analysis_counter; + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Pick up nodes*/ + if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17472) @@ -53,7 +53,4 @@ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); #endif - void InputUpdateFromConstant(IssmDouble constant, int name); - void InputUpdateFromConstant(int constant, int name); - void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ @@ -79,6 +76,4 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); int GetNodeIndex(Node* node); - void GetNodesSidList(int* sidlist); - void GetNodesLidList(int* lidlist); int GetNumberOfNodes(void); int GetNumberOfNodesPressure(void); @@ -102,8 +97,6 @@ void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); void GetVectorFromInputs(Vector* vector, int name_enum); - void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); - void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputDuplicate(int original_enum,int new_enum); @@ -215,5 +208,4 @@ void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); Node* GetNode(int node_number); - void InputChangeName(int enum_type,int enum_type_old); void InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type); void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");}; @@ -232,7 +224,4 @@ Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order); - ElementVector* NewElementVector(int approximation_enum); - ElementMatrix* NewElementMatrix(int approximation_enum); - ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss); void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp (revision 17472) +++ /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp (revision 17472) @@ -0,0 +1,419 @@ +/*!\file TetraInput.c + * \brief: implementation of the TetraInput object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../classes.h" +#include "../../shared/shared.h" + +/*TetraInput constructors and destructor*/ +/*FUNCTION TetraInput::TetraInput(){{{*/ +TetraInput::TetraInput(){ + values = NULL; +} +/*}}}*/ +/*FUNCTION TetraInput::TetraInput(int in_enum_type,IssmDouble* invalues,element_type_in){{{*/ +TetraInput::TetraInput(int in_enum_type,IssmDouble* in_values,int element_type_in) + :TetraRef(1) +{ + + /*Set TetraRef*/ + this->SetElementType(element_type_in,0); + this->element_type=element_type_in; + + /*Set Enum*/ + enum_type=in_enum_type; + + /*Set values*/ + this->values=xNew(this->NumberofNodes()); + for(int i=0;iNumberofNodes();i++) values[i]=in_values[i]; +} +/*}}}*/ +/*FUNCTION TetraInput::~TetraInput(){{{*/ +TetraInput::~TetraInput(){ + xDelete(this->values); +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION TetraInput::Echo {{{*/ +void TetraInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION TetraInput::DeepEcho{{{*/ +void TetraInput::DeepEcho(void){ + + _printf_(setw(15)<<" TetraInput "<enum_type)<<" ["); + for(int i=0;iNumberofNodes();i++) _printf_(" "<values[i]); + _printf_("] ("<element_type)<<")\n"); +} +/*}}}*/ +/*FUNCTION TetraInput::Id{{{*/ +int TetraInput::Id(void){ return -1; } +/*}}}*/ +/*FUNCTION TetraInput::ObjectEnum{{{*/ +int TetraInput::ObjectEnum(void){ + + return TetraInputEnum; + +} +/*}}}*/ +/*FUNCTION TetraInput::copy{{{*/ +Object* TetraInput::copy() { + + return new TetraInput(this->enum_type,this->values,this->element_type); + +} +/*}}}*/ + +/*TetraInput management*/ +/*FUNCTION TetraInput::InstanceEnum{{{*/ +int TetraInput::InstanceEnum(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION TetraInput::GetResultInterpolation{{{*/ +int TetraInput::GetResultInterpolation(void){ + + return P1Enum; + +} +/*}}}*/ +/*FUNCTION TetraInput::GetResultNumberOfNodes{{{*/ +int TetraInput::GetResultNumberOfNodes(void){ + + return this->NumberofNodes(); + +} +/*}}}*/ +/*FUNCTION TetraInput::ResultToPatch{{{*/ +void TetraInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){ + + int numnodes = this->NumberofNodes(); + + /*Some checks*/ + _assert_(values); + _assert_(numnodes==nodesperelement); + + /*Fill in arrays*/ + for(int i=0;ivalues[i]; +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION TetraInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){{{*/ +void TetraInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){ + + /*Call TetraRef function*/ + _assert_(gauss->Enum()==GaussTetraEnum); + TetraRef::GetInputValue(pvalue,&values[0],(GaussTetra*)gauss); + +} +/*}}}*/ +/*FUNCTION TetraInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){{{*/ +void TetraInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){ + + /*Call TetraRef function*/ + _assert_(gauss->Enum()==GaussTetraEnum); + TetraRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTetra*)gauss); +} +/*}}}*/ +/*FUNCTION TetraInput::ChangeEnum{{{*/ +void TetraInput::ChangeEnum(int newenumtype){ + this->enum_type=newenumtype; +} +/*}}}*/ +/*FUNCTION TetraInput::GetInputAverage{{{*/ +void TetraInput::GetInputAverage(IssmDouble* pvalue){ + + int numnodes = this->NumberofNodes(); + IssmDouble numnodesd = reCast(numnodes); + IssmDouble value = 0.; + + for(int i=0;i(1); + times=xNew(1); + + this->GetInputAverage(&outvalues[0]); + times[0]=0.; /*we don't have a time*/ + + *pvalues=outvalues; + *ptimes=times; + *pnumtimes=numtimes; +} +/*}}}*/ +/*FUNCTION TetraInput::GetInputUpToCurrentTimeAverages{{{*/ +void TetraInput::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){ + + IssmDouble* outvalues=NULL; + IssmDouble* times=NULL; + int numtimes; + + /*this is not a transient forcing, so we only have 1 value, steady state: */ + numtimes=1; + outvalues=xNew(1); + times=xNew(1); + + this->GetInputAverage(&outvalues[0]); + times[0]=currenttime; /*we don't have a time*/ + + *pvalues=outvalues; + *ptimes=times; + *pnumtimes=numtimes; +} +/*}}}*/ + +/*Intermediary*/ +/*FUNCTION TetraInput::SquareMin{{{*/ +void TetraInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){ + + int numnodes=this->NumberofNodes(); + IssmDouble squaremin; + + /*Now, figure out minimum of valuescopy: */ + squaremin=pow(this->values[0],2); + for(int i=1;ivalues[i],2)values[i],2); + } + /*Assign output pointers:*/ + *psquaremin=squaremin; +} +/*}}}*/ +/*FUNCTION TetraInput::ContrainMin{{{*/ +void TetraInput::ConstrainMin(IssmDouble minimum){ + + int numnodes = this->NumberofNodes(); + for(int i=0;iNumberofNodes(); + + for(int i=0;inorm) norm=fabs(values[i]); + return norm; +} +/*}}}*/ +/*FUNCTION TetraInput::Max{{{*/ +IssmDouble TetraInput::Max(void){ + + int numnodes=this->NumberofNodes(); + IssmDouble max=values[0]; + + for(int i=1;imax) max=values[i]; + } + return max; +} +/*}}}*/ +/*FUNCTION TetraInput::MaxAbs{{{*/ +IssmDouble TetraInput::MaxAbs(void){ + + int numnodes=this->NumberofNodes(); + IssmDouble max=fabs(values[0]); + + for(int i=1;imax) max=fabs(values[i]); + } + return max; +} +/*}}}*/ +/*FUNCTION TetraInput::Min{{{*/ +IssmDouble TetraInput::Min(void){ + + const int numnodes=this->NumberofNodes(); + IssmDouble min=values[0]; + + for(int i=1;iNumberofNodes(); + IssmDouble min=fabs(values[0]); + + for(int i=1;iNumberofNodes(); + for(int i=0;iNumberofNodes(); + for(int i=0;iNumberofNodes(); + TetraInput* xtriainput=NULL; + + /*xinput is of the same type, so cast it: */ + if(xinput->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + xtriainput=(TetraInput*)xinput; + if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + + /*Carry out the AXPY operation depending on type:*/ + for(int i=0;ivalues[i]=this->values[i]+scalar*xtriainput->values[i]; + +} +/*}}}*/ +/*FUNCTION TetraInput::Constrain{{{*/ +void TetraInput::Constrain(IssmDouble cm_min, IssmDouble cm_max){ + + int i; + const int numnodes=this->NumberofNodes(); + + if(!xIsNan(cm_min)) for(i=0;ivalues[i]values[i]=cm_min; + if(!xIsNan(cm_max)) for(i=0;ivalues[i]>cm_max)this->values[i]=cm_max; + +} +/*}}}*/ +/*FUNCTION TetraInput::GetVectorFromInputs{{{*/ +void TetraInput::GetVectorFromInputs(Vector* vector,int* doflist){ + + const int numnodes=this->NumberofNodes(); + vector->SetValues(numnodes,doflist,this->values,INS_VAL); + +} /*}}}*/ +/*FUNCTION TetraInput::PointwiseMin{{{*/ +Input* TetraInput::PointwiseMin(Input* inputB){ + + /*Ouput*/ + TetraInput* outinput=NULL; + + /*Intermediaries*/ + int i; + TetraInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *minvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(TetraInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*Create point wise min*/ + for(i=0;ivalues[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i]; + else minvalues[i]=this->values[i]; + } + + /*Create new Tetra vertex input (copy of current input)*/ + outinput=new TetraInput(this->enum_type,&minvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(minvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION TetraInput::PointwiseMax{{{*/ +Input* TetraInput::PointwiseMax(Input* inputB){ + + /*Ouput*/ + TetraInput* outinput=NULL; + + /*Intermediaries*/ + int i; + TetraInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *maxvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(TetraInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*Create point wise max*/ + for(i=0;ivalues[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i]; + else maxvalues[i]=this->values[i]; + } + + /*Create new Tetra vertex input (copy of current input)*/ + outinput=new TetraInput(this->enum_type,&maxvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(maxvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION TetraInput::PointwiseDivide{{{*/ +Input* TetraInput::PointwiseDivide(Input* inputB){ + + /*Ouput*/ + TetraInput* outinput=NULL; + + /*Intermediaries*/ + TetraInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(TetraInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*Allocate intermediary*/ + IssmDouble* AdotBvalues=xNew(numnodes); + + /*Create point wise division*/ + for(int i=0;ivalues[i]!=0); + AdotBvalues[i]=this->values[i]/xinputB->values[i]; + } + + /*Create new Tetra vertex input (copy of current input)*/ + outinput=new TetraInput(this->enum_type,AdotBvalues,this->element_type); + + /*Return output pointer*/ + xDelete(AdotBvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION TetraInput::Configure{{{*/ +void TetraInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h (revision 17472) +++ /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h (revision 17472) @@ -0,0 +1,75 @@ +/*! \file TetraInput.h + * \brief: header file for TetraInput object + */ + +#ifndef _TETRAINPUT_H_ +#define _TETRAINPUT_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Input.h" +#include "../Elements/TetraRef.h" +class Gauss; +class Gauss; +/*}}}*/ + +class TetraInput: public Input,public TetraRef{ + + public: + int enum_type; + IssmDouble* values; + + /*TetraInput constructors, destructors*/ + TetraInput(); + TetraInput(int enum_type,IssmDouble* values,int element_type_in); + ~TetraInput(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho(); + int Id(); + int ObjectEnum(); + Object *copy(); + + /*TetraInput management:*/ + int InstanceEnum(); + Input* SpawnTriaInput(int location){_error_("not supported yet");}; + Input* SpawnSegInput(int index1,int index2){_error_("not supported yet");}; + Input* PointwiseDivide(Input* inputB); + Input* PointwiseMin(Input* inputB); + Input* PointwiseMax(Input* inputB); + int GetResultInterpolation(void); + int GetResultNumberOfNodes(void); + void ResultToPatch(IssmDouble* values,int nodesperelement,int sid); + void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; + void Configure(Parameters* parameters); + + /*numerics*/ + void GetInputValue(bool* pvalue){_error_("not implemented yet");} + void GetInputValue(int* pvalue){_error_("not implemented yet");} + void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");} + void GetInputValue(IssmDouble* pvalue,Gauss* gauss); + void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int index){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss); + void GetInputAverage(IssmDouble* pvalue); + void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes); + void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime); + void ChangeEnum(int newenumtype); + void SquareMin(IssmDouble* psquaremin,Parameters* parameters); + void ConstrainMin(IssmDouble minimum); + void Set(IssmDouble setvalue); + void Scale(IssmDouble scale_factor); + void AXPY(Input* xinput,IssmDouble scalar); + void Constrain(IssmDouble cm_min, IssmDouble cm_max); + IssmDouble InfinityNorm(void); + IssmDouble Max(void); + IssmDouble MaxAbs(void); + IssmDouble Min(void); + IssmDouble MinAbs(void); + void Extrude(void){_error_("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");}; + void GetVectorFromInputs(Vector* vector,int* doflist); + +}; +#endif /* _TETRAINPUT_H */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 17471) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 17472) @@ -61,4 +61,5 @@ #include "./Inputs/DoubleInput.h" #include "./Inputs/IntInput.h" +#include "./Inputs/TetraInput.h" #include "./Inputs/PentaInput.h" #include "./Inputs/TriaInput.h" Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp (revision 17472) @@ -3,4 +3,5 @@ */ +#include #include "./GaussTetra.h" #include "../../shared/io/Print/Print.h" @@ -32,5 +33,42 @@ /*FUNCTION GaussTetra::GaussTetra(int order) {{{*/ GaussTetra::GaussTetra(int order){ - _error_("not implemented yet"); + /*Get gauss points*/ + GaussLegendreTetra(&numgauss,&coords1,&coords2,&coords3,&coords4,&weights,order); + for(int i=0;iweights[i]=this->weights[i]*sqrt(2.)/6.; //FIXME: double check that + + /*Initialize static fields as undefinite*/ + weight=UNDEF; + coord1=UNDEF; + coord2=UNDEF; + coord3=UNDEF; +} +/*}}}*/ +/*FUNCTION GaussTetra::GaussTetra(int index1,int index2,int index3,int order) {{{*/ +GaussTetra::GaussTetra(int index1,int index2,int index3,int order){ + + /*Basal Tria*/ + if(index1==0 && index2==1 && index3==2){ + GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order); + coords4=xNew(numgauss); + for(int i=0;i(numgauss); + for(int i=0;i(numgauss); + for(int i=0;i(numgauss); + for(int i=0;i(numberofelements); - npart=xNew(numberofnodes); - index=xNew(elements_width*numberofelements); - for (i=0;i1){ -#ifdef _HAVE_METIS_ - METIS_PartMeshNodalPatch(&numberofelements,&numberofnodes, index, &etype, &numflag, &num_procs, &edgecut, epart, npart); -#else - _error_("metis has not beed installed. Cannot run with more than 1 cpu"); -#endif - } - else if (num_procs==1){ - /*METIS does not know how to deal with one cpu only!*/ - for (i=0;i(numberofelements); npart=xNew(numberofnodes); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 17472) @@ -37,4 +37,9 @@ } break; + case Mesh3DtetrasEnum: + for(i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]) elements->AddObject(new Tetra(i+1,i,i,iomodel,nummodels)); + } + break; case Mesh3DEnum: iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum); @@ -60,5 +65,5 @@ if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBbarEnum,QmuMaterialsRheologyBEnum); break; - case Mesh3DEnum: + case Mesh3DEnum: case Mesh3DtetrasEnum: if(dakota_analysis) elements->InputDuplicate(MaterialsRheologyBEnum,QmuMaterialsRheologyBEnum); break; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 17472) @@ -39,4 +39,5 @@ case Mesh2DhorizontalEnum: elementswidth=3; break; case Mesh2DverticalEnum: elementswidth=3; break; + case Mesh3DtetrasEnum: elementswidth=4; break; case Mesh3DEnum: elementswidth=6; break; default: _error_("mesh not supported yet"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 17472) @@ -60,4 +60,10 @@ case Mesh2DverticalEnum: elements_width=3; + numberofelements2d = 0; + numberofvertices2d = 0; + numlayers = 0; + break; + case Mesh3DtetrasEnum: + elements_width=4; //penta elements numberofelements2d = 0; numberofvertices2d = 0; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17471) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17472) @@ -215,4 +215,5 @@ Mesh2DverticalEnum, Mesh3DEnum, + Mesh3DtetrasEnum, MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script) MasstransportHydrostaticAdjustmentEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17472) @@ -223,4 +223,5 @@ case Mesh2DverticalEnum : return "Mesh2Dvertical"; case Mesh3DEnum : return "Mesh3D"; + case Mesh3DtetrasEnum : return "Mesh3Dtetras"; case MiscellaneousNameEnum : return "MiscellaneousName"; case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17471) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17472) @@ -226,4 +226,5 @@ else if (strcmp(name,"Mesh2Dvertical")==0) return Mesh2DverticalEnum; else if (strcmp(name,"Mesh3D")==0) return Mesh3DEnum; + else if (strcmp(name,"Mesh3Dtetras")==0) return Mesh3DtetrasEnum; else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum; else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum; - else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; + if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; + else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Nodes")==0) return NodesEnum; - else if (strcmp(name,"Contours")==0) return ContoursEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Parameters")==0) return ParametersEnum; + if (strcmp(name,"Contours")==0) return ContoursEnum; + else if (strcmp(name,"Parameters")==0) return ParametersEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"Results")==0) return ResultsEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"Vz")==0) return VzEnum; else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; - else if (strcmp(name,"VzHO")==0) return VzHOEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"VzPicard")==0) return VzPicardEnum; + if (strcmp(name,"VzHO")==0) return VzHOEnum; + else if (strcmp(name,"VzPicard")==0) return VzPicardEnum; else if (strcmp(name,"VzFS")==0) return VzFSEnum; else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum; - else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; + if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; + else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;