Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16889) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16890) @@ -27,5 +27,12 @@ /*Finite Element Analysis*/ ElementMatrix* AdjointBalancethicknessAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + BalancethicknessAnalysis* analysis = new BalancethicknessAnalysis(); + ElementMatrix* Ke = analysis->CreateKMatrix(element); + delete analysis; + + /*Transpose and return Ke*/ + Ke->Transpose(); + return Ke; }/*}}}*/ ElementVector* AdjointBalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16889) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16890) @@ -38,9 +38,163 @@ /*Finite Element Analysis*/ ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixVolume(element); + ElementMatrix* Ke2=CreateKMatrixSurface(element); + ElementMatrix* Ke3=CreateKMatrixBed(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble Jdet,D; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* B = xNew(numnodes); + IssmDouble* Bprime = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(Bprime,gauss); + GetB(B,element,xyz_list,gauss); + D=gauss->weight*Jdet; + + TripleMultiply(B,1,numnodes,1, + &D,1,1,0, + Bprime,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; +} +/*}}}*/ +ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ + + if(!element->IsOnSurface()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet,D,normal[2]; + IssmDouble *xyz_list_top = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*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[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->IsOnBed()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet,D,normal[2]; + IssmDouble *xyz_list_base = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*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[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,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.*/ + + /*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 ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 16889) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 16890) @@ -22,5 +22,9 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixVolume(Element* element); + ElementMatrix* CreateKMatrixSurface(Element* element); + ElementMatrix* CreateKMatrixBed(Element* element); ElementVector* CreatePVector(Element* element); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16889) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16890) @@ -66,4 +66,8 @@ case Mesh2DhorizontalEnum: basalelement = element; + break; + case Mesh2DverticalEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); break; case Mesh3DEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16889) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16890) @@ -1213,4 +1213,22 @@ }/*}}}*/ +/*FUNCTION Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/ +void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){ + + int indices[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Allocate Output*/ + IssmDouble* xyz_list_edge = xNew(2*3); + this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); + for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; + + /*Assign output pointer*/ + *pxyz_list = xyz_list_edge; + +}/*}}}*/ /*FUNCTION Tria::NormalSection{{{*/ void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){ @@ -2086,4 +2104,12 @@ } /*}}}*/ +/*FUNCTION Tria::JacobianDeterminantTop{{{*/ +void Tria::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetSegmentJacobianDeterminant(pJdet,xyz_list_top,(GaussTria*)gauss); + +} +/*}}}*/ /*FUNCTION Tria::HasEdgeOnBed {{{*/ bool Tria::HasEdgeOnBed(){ @@ -2310,4 +2336,12 @@ int indices[2]; this->EdgeOnBedIndices(&indices[0],&indices[1]); + return new GaussTria(indices[0],indices[1],order); +} +/*}}}*/ +/*FUNCTION Tria::NewGaussTop(int order){{{*/ +Gauss* Tria::NewGaussTop(int order){ + + int indices[2]; + this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); return new GaussTria(indices[0],indices[1],order); } @@ -2412,4 +2446,20 @@ /*FUNCTION Tria::NormalBase {{{*/ void Tria::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ + + /*Build unit outward pointing vector*/ + IssmDouble vector[2]; + IssmDouble norm; + + vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; + vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; + + norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]); + + bed_normal[0]= + vector[1]/norm; + bed_normal[1]= - vector[0]/norm; +} +/*}}}*/ +/*FUNCTION Tria::NormalTop {{{*/ +void Tria::NormalTop(IssmDouble* bed_normal,IssmDouble* xyz_list){ /*Build unit outward pointing vector*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16889) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16890) @@ -121,5 +121,5 @@ void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); - void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; + 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,int object_enum=MeshElementsEnum); @@ -270,5 +270,5 @@ IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented");}; void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); - void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; + void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); IssmDouble GetMaterialParameter(int enum_in); @@ -294,5 +294,5 @@ void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); - void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; Gauss* NewGauss(void); @@ -303,5 +303,5 @@ Gauss* NewGaussBase(int order); Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; - Gauss* NewGaussTop(int order){_error_("not implemented yet");}; + Gauss* NewGaussTop(int order); ElementVector* NewElementVector(int approximation_enum); ElementMatrix* NewElementMatrix(int approximation_enum);