Index: /issm/trunk-jpl/src/c/analyses/Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 18929) @@ -26,10 +26,10 @@ /*Model processing*/ - virtual int DofsPerNode(int** doflist,int domaintype,int approximation)=0; - virtual void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum)=0; - virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0; - virtual void CreateNodes(Nodes* nodes,IoModel* iomodel)=0; virtual void CreateConstraints(Constraints* constraints,IoModel* iomodel)=0; virtual void CreateLoads(Loads* loads, IoModel* iomodel)=0; + virtual void CreateNodes(Nodes* nodes,IoModel* iomodel)=0; + virtual int DofsPerNode(int** doflist,int domaintype,int approximation)=0; + virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0; + virtual void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum)=0; /*Finite element Analysis*/ @@ -39,8 +39,8 @@ virtual ElementMatrix* CreateKMatrix(Element* element)=0; virtual ElementVector* CreatePVector(Element* element)=0; - virtual void GetSolutionFromInputs(Vector* solution,Element* element)=0; - virtual void GradientJ(Vector* gradient,Element* element,int control_type,int control_index)=0; - virtual void InputUpdateFromSolution(IssmDouble* solution,Element* element)=0; - virtual void UpdateConstraints(FemModel* femmodel)=0; + virtual void GetSolutionFromInputs(Vector* solution,Element* element)=0; + virtual void GradientJ(Vector* gradient,Element* element,int control_type,int control_index)=0; + virtual void InputUpdateFromSolution(IssmDouble* solution,Element* element)=0; + virtual void UpdateConstraints(FemModel* femmodel)=0; }; #endif Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp (revision 18929) @@ -6,8 +6,15 @@ /*Model processing*/ +void DepthAverageAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void DepthAverageAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void DepthAverageAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(nodes,iomodel,DepthAverageAnalysisEnum,P1Enum); + +}/*}}}*/ int DepthAverageAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void DepthAverageAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void DepthAverageAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -26,12 +33,5 @@ } }/*}}}*/ -void DepthAverageAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - ::CreateNodes(nodes,iomodel,DepthAverageAnalysisEnum,P1Enum); - -}/*}}}*/ -void DepthAverageAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ -}/*}}}*/ -void DepthAverageAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +void DepthAverageAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -137,5 +137,5 @@ return pe; }/*}}}*/ -void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; where hi is the interpolation function for node i.*/ @@ -157,11 +157,11 @@ } /*}}}*/ -void DepthAverageAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void DepthAverageAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void DepthAverageAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void DepthAverageAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void DepthAverageAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void DepthAverageAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int inputenum; @@ -169,5 +169,5 @@ element->InputUpdateFromSolutionOneDof(solution,inputenum); }/*}}}*/ -void DepthAverageAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void DepthAverageAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,9 +26,9 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 18929) @@ -6,11 +6,22 @@ #include "../solutionsequences/solutionsequences.h" -int ExtrapolationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ +void ExtrapolationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + // do nothing for now + return; +} +/*}}}*/ +void ExtrapolationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + // do nothing for now + return; +}/*}}}*/ +void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + int finiteelement=P1Enum; + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement); + iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); +} +/*}}}*/ +int ExtrapolationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -} -/*}}}*/ -void ExtrapolationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ - //do nothing for now - return; } /*}}}*/ @@ -36,23 +47,12 @@ } /*}}}*/ -void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - int finiteelement=P1Enum; - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); - ::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement); - iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); -} -/*}}}*/ -void ExtrapolationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - // do nothing for now +void ExtrapolationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + //do nothing for now return; } /*}}}*/ -void ExtrapolationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - // do nothing for now - return; -}/*}}}*/ /*Finite element Analysis*/ -void ExtrapolationAnalysis::Core(FemModel* femmodel){/*{{{*/ +void ExtrapolationAnalysis::Core(FemModel* femmodel){/*{{{*/ /* Intermediaries */ @@ -237,11 +237,66 @@ return pe; }/*}}}*/ -void ExtrapolationAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For node i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ N ] + * [ N ] + * where N is the finiteelement function for node i. + * + * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void 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); + +}/*}}}*/ +void ExtrapolationAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void ExtrapolationAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void ExtrapolationAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int domaintype, extrapolationvariable; @@ -258,60 +313,5 @@ } }/*}}}*/ -void ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. - * For node i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ N ] - * [ N ] - * where N is the finiteelement function for node i. - * - * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - element->NodalFunctions(basis,gauss); - - /*Build B: */ - for(int i=0;i(basis); -}/*}}}*/ -void 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); - -}/*}}}*/ -void ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/ +void ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/ int numnodes=element->GetNumberOfNodes(); @@ -347,5 +347,5 @@ delete gauss; }/*}}}*/ -void ExtrapolationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void ExtrapolationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ for(int i=0;ielements->Size();i++){ Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,11 +26,11 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); - void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void SetConstraintsOnIce(Element* element); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void SetConstraintsOnIce(Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 18929) @@ -6,8 +6,15 @@ /*Model processing*/ +void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(nodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum); + +}/*}}}*/ int ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void ExtrudeFromBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -26,12 +33,5 @@ } }/*}}}*/ -void ExtrudeFromBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - ::CreateNodes(nodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum); - -}/*}}}*/ -void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ -}/*}}}*/ -void ExtrudeFromBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -61,4 +61,94 @@ return Ke; }/*}}}*/ +ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/ + + if(!element->IsOnBase()) return NULL; + + /*Intermediaries */ + int dim; + IssmDouble Jdet,D,normal[3]; + IssmDouble *xyz_list_base = NULL; + + /*Get dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->NormalBase(&normal[0],xyz_list_base); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctions(basis,gauss); + D = - gauss->weight*Jdet*normal[dim-1]; + + TripleMultiply(basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list_base); + xDelete(basis); + delete gauss; + return Ke; +} +/*}}}*/ +ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ + + if(!element->IsOnSurface()) return NULL; + + /*Intermediaries */ + int dim; + IssmDouble Jdet,D,normal[3]; + IssmDouble *xyz_list_top = NULL; + + /*Get dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesTop(&xyz_list_top); + element->NormalTop(&normal[0],xyz_list_top); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussTop(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss); + element->NodalFunctions(basis,gauss); + D = - gauss->weight*Jdet*normal[dim-1]; + + TripleMultiply(basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list_top); + xDelete(basis); + delete gauss; + return Ke; +} +/*}}}*/ ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/ @@ -106,98 +196,8 @@ } /*}}}*/ -ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ - - if(!element->IsOnSurface()) return NULL; - - /*Intermediaries */ - int dim; - IssmDouble Jdet,D,normal[3]; - IssmDouble *xyz_list_top = NULL; - - /*Get dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesTop(&xyz_list_top); - element->NormalTop(&normal[0],xyz_list_top); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussTop(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss); - element->NodalFunctions(basis,gauss); - D = - gauss->weight*Jdet*normal[dim-1]; - - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(xyz_list_top); - xDelete(basis); - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/ - - if(!element->IsOnBase()) return NULL; - - /*Intermediaries */ - int dim; - IssmDouble Jdet,D,normal[3]; - IssmDouble *xyz_list_base = NULL; - - /*Get dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->NormalBase(&normal[0],xyz_list_base); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctions(basis,gauss); - D = - gauss->weight*Jdet*normal[dim-1]; - - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(xyz_list_base); - xDelete(basis); - delete gauss; - return Ke; -} -/*}}}*/ ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/ return NULL; }/*}}}*/ -void ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; where hi is the interpolation function for node i.*/ @@ -219,11 +219,11 @@ } /*}}}*/ -void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void ExtrudeFromBaseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void ExtrudeFromBaseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void ExtrudeFromBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void ExtrudeFromBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int inputenum; @@ -231,5 +231,5 @@ element->InputUpdateFromSolutionOneDof(solution,inputenum); }/*}}}*/ -void ExtrudeFromBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void ExtrudeFromBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -25,13 +25,13 @@ ElementMatrix* CreateJacobianMatrix(Element* element); ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixBed(Element* element); + ElementMatrix* CreateKMatrixSurface(Element* element); ElementMatrix* CreateKMatrixVolume(Element* element); - ElementMatrix* CreateKMatrixSurface(Element* element); - ElementMatrix* CreateKMatrixBed(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp (revision 18929) @@ -6,8 +6,15 @@ /*Model processing*/ +void ExtrudeFromTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(nodes,iomodel,ExtrudeFromTopAnalysisEnum,P1Enum); + +}/*}}}*/ int ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void ExtrudeFromTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void ExtrudeFromTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -26,12 +33,5 @@ } }/*}}}*/ -void ExtrudeFromTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - ::CreateNodes(nodes,iomodel,ExtrudeFromTopAnalysisEnum,P1Enum); - -}/*}}}*/ -void ExtrudeFromTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ -}/*}}}*/ -void ExtrudeFromTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +void ExtrudeFromTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -199,5 +199,5 @@ return NULL; }/*}}}*/ -void ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; where hi is the interpolation function for node i.*/ @@ -219,11 +219,11 @@ } /*}}}*/ -void ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void ExtrudeFromTopAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void ExtrudeFromTopAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void ExtrudeFromTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void ExtrudeFromTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int inputenum; @@ -231,5 +231,5 @@ element->InputUpdateFromSolutionOneDof(solution,inputenum); }/*}}}*/ -void ExtrudeFromTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void ExtrudeFromTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -25,13 +25,13 @@ ElementMatrix* CreateJacobianMatrix(Element* element); ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixBed(Element* element); + ElementMatrix* CreateKMatrixSurface(Element* element); ElementMatrix* CreateKMatrixVolume(Element* element); - ElementMatrix* CreateKMatrixSurface(Element* element); - ElementMatrix* CreateKMatrixBed(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 18929) @@ -6,8 +6,55 @@ /*Model processing*/ +void FreeSurfaceBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void FreeSurfaceBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonbase=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + if(!(reCast(nodeonbase[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbase[reCast(vertex_pairing[2*i+1])-1]))) continue; + } + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceBaseAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); +}/*}}}*/ +void FreeSurfaceBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(nodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ int FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void FreeSurfaceBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -43,52 +90,5 @@ } }/*}}}*/ -void FreeSurfaceBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); - ::CreateNodes(nodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum); - iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); -}/*}}}*/ -void FreeSurfaceBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ -}/*}}}*/ -void FreeSurfaceBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - - /*Intermediaries*/ - int penpair_ids[2]; - int count=0; - int numvertex_pairing; - - /*Create Penpair for vertex_pairing: */ - IssmDouble *vertex_pairing=NULL; - IssmDouble *nodeonbase=NULL; - iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); - for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ - - /*In debugging mode, check that the second node is in the same cpu*/ - _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); - - /*Skip if one of the two is not on the bed*/ - if(iomodel->domaintype!=Domain2DhorizontalEnum){ - if(!(reCast(nodeonbase[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbase[reCast(vertex_pairing[2*i+1])-1]))) continue; - } - - /*Get node ids*/ - penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); - penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); - - /*Create Load*/ - loads->AddObject(new Penpair( - iomodel->loadcounter+count+1, - &penpair_ids[0], - FreeSurfaceBaseAnalysisEnum)); - count++; - } - } - - /*free ressources: */ - iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); - iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); +void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -304,5 +304,5 @@ }/*}}}*/ -void FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,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 @@ -332,5 +332,5 @@ xDelete(basis); }/*}}}*/ -void FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,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 @@ -361,14 +361,14 @@ }/*}}}*/ -void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void FreeSurfaceBaseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void FreeSurfaceBaseAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void FreeSurfaceBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void FreeSurfaceBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ element->InputUpdateFromSolutionOneDof(solution,BaseEnum); }/*}}}*/ -void FreeSurfaceBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void FreeSurfaceBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Intermediary*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,10 +26,10 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 18929) @@ -6,8 +6,55 @@ /*Model processing*/ +void FreeSurfaceTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void FreeSurfaceTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonsurface=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + if(!(reCast(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; + } + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceTopAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum); +}/*}}}*/ +void FreeSurfaceTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(nodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ int FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void FreeSurfaceTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void FreeSurfaceTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -51,52 +98,5 @@ } }/*}}}*/ -void FreeSurfaceTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); - ::CreateNodes(nodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum); - iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum); -}/*}}}*/ -void FreeSurfaceTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ -}/*}}}*/ -void FreeSurfaceTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - - /*Intermediaries*/ - int penpair_ids[2]; - int count=0; - int numvertex_pairing; - - /*Create Penpair for vertex_pairing: */ - IssmDouble *vertex_pairing=NULL; - IssmDouble *nodeonsurface=NULL; - iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); - for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ - - /*In debugging mode, check that the second node is in the same cpu*/ - _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); - - /*Skip if one of the two is not on the bed*/ - if(iomodel->domaintype!=Domain2DhorizontalEnum){ - if(!(reCast(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; - } - - /*Get node ids*/ - penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); - penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); - - /*Create Load*/ - loads->AddObject(new Penpair( - iomodel->loadcounter+count+1, - &penpair_ids[0], - FreeSurfaceTopAnalysisEnum)); - count++; - } - } - - /*free ressources: */ - iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); - iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum); +void FreeSurfaceTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -307,5 +307,5 @@ }/*}}}*/ -void FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,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 @@ -335,5 +335,5 @@ xDelete(basis); }/*}}}*/ -void FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,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 @@ -364,15 +364,15 @@ }/*}}}*/ -void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void FreeSurfaceTopAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void FreeSurfaceTopAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void FreeSurfaceTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void FreeSurfaceTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum); }/*}}}*/ -void FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,10 +26,10 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp (revision 18929) @@ -6,8 +6,15 @@ /*Model processing*/ +void GiaAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + /*No constraints*/ +}/*}}}*/ +void GiaAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ +void GiaAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + ::CreateNodes(nodes,iomodel,GiaAnalysisEnum,P1Enum); +}/*}}}*/ int GiaAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; -}/*}}}*/ -void GiaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ void GiaAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -27,12 +34,5 @@ iomodel->FetchDataToInput(elements,GiaLithosphereThicknessEnum); }/*}}}*/ -void GiaAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - ::CreateNodes(nodes,iomodel,GiaAnalysisEnum,P1Enum); -}/*}}}*/ -void GiaAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - /*No constraints*/ -}/*}}}*/ -void GiaAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - /*No loads*/ +void GiaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ }/*}}}*/ @@ -54,14 +54,14 @@ _error_("not implemented yet"); }/*}}}*/ -void GiaAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void GiaAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void GiaAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void GiaAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void GiaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void GiaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ _error_("not implemented yet"); }/*}}}*/ -void GiaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void GiaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h (revision 18928) +++ /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h (revision 18929) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,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); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,8 +26,8 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); }; #endif