Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 10366) +++ /issm/trunk/src/c/Makefile.am (revision 10367) @@ -184,4 +184,5 @@ ./shared/Numerics/cross.cpp\ ./shared/Numerics/extrema.cpp\ + ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ ./shared/Numerics/UnitConversion.cpp\ ./shared/Numerics/PetscOptionsFromAnalysis.cpp\ @@ -198,4 +199,5 @@ ./shared/Elements/GetGlobalDofList.cpp\ ./shared/Elements/GetNumberOfDofs.cpp\ + ./shared/Elements/CoordinateSystemTransform.cpp\ ./shared/String/sharedstring.h\ ./toolkits/petsc\ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 10366) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 10367) @@ -2830,4 +2830,5 @@ ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction(); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + TransformStiffnessMatrixCoord(Ke,2); /*clean-up and return*/ @@ -3031,4 +3032,7 @@ } + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,2); + /*Clean up and return*/ delete gauss; @@ -3188,6 +3192,10 @@ /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;inrows*Ke->ncols*sizeof(double)); + for(i=0;inrows;i++) for(j=0;jncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes,numnodes,dim); + + /*Transform matrix: T*Ke*T^t */ + TripleMultiply(transform,numnodes*dim,numnodes*dim,1, + values,Ke->nrows,Ke->ncols,0, + transform,numnodes*dim,numnodes*dim,0, + &Ke->values[0],0); + + /*Free Matrix*/ + xfree((void**)&transform); + xfree((void**)&values); +} +/*}}}*/ +/*FUNCTION Tria::TransformLoadVectorCoord{{{1*/ +void Tria::TransformLoadVectorCoord(ElementVector* pe,int dim){ + + int i,j; + int numnodes = NUMVERTICES; + double *transform = NULL; + double *values = NULL; + + /*Copy current load vector*/ + values=(double*)xmalloc(pe->nrows*sizeof(double)); + for(i=0;inrows;i++) values[i]=pe->values[i]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes,numnodes,dim); + + /*Transform matrix: T*pe */ + MatrixMultiply(transform,numnodes*dim,numnodes*dim,1, + values,pe->nrows,1,0, + &pe->values[0],0); + + /*Free Matrix*/ + xfree((void**)&transform); + xfree((void**)&values); +} +/*}}}*/ +/*FUNCTION Tria::TransformSolutionCoord{{{1*/ +void Tria::TransformSolutionCoord(double* vel,int dim,int node_index){ + + int i,j; + double *transform = NULL; + double *values = NULL; + + /*Copy current solution vector*/ + values=(double*)xmalloc(dim*sizeof(double)); + for(i=0;inrows*sizeof(double)); + for(int i=0;inrows;i++) values[i]=pe->values[i]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes,numnodes,dim); + + /*Transform matrix: T*pe */ + MatrixMultiply(transform,numnodes*dim,numnodes*dim,1, + values,pe->nrows,1,0, + &pe->values[0],0); + + /*Free Matrix*/ + xfree((void**)&transform); + xfree((void**)&values); } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 10366) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 10367) @@ -80,7 +80,8 @@ /*}}}*/ /*Load management: {{{1*/ - void GetDofList(int** pdoflist,int approximation_enum,int setenum); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetSegmentNormal(double* normal,double xyz_list[2][3]); void GetQuadNormal(double* normal,double xyz_list[4][3]); + void TransformLoadVectorCoord(ElementVector* pe,int dim); #ifdef _HAVE_CONTROL_ ElementVector* CreatePVectorAdjointHoriz(void); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 10366) +++ /issm/trunk/src/c/objects/Node.cpp (revision 10367) @@ -43,5 +43,8 @@ this->sid=node_sid; this->analysis_type=analysis_type; - for(k=0;k<3;k++) for(l=0;l<3;l++) this->referential[k][l]=1.0; + + /*Initialize coord_system: Identity matrix by default*/ + for(k=0;k<3;k++) for(l=0;l<3;l++) this->coord_system[k][l]=0.0; + for(k=0;k<3;k++) this->coord_system[k][k]=1.0; /*indexing:*/ @@ -54,10 +57,16 @@ //intialize inputs, and add as many inputs per element as requested: this->inputs=new Inputs(); - if (iomodel->Data(MeshVertexonbedEnum)) this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index])); - if (iomodel->Data(MeshVertexonsurfaceEnum)) this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index])); - if (iomodel->Data(MaskVertexonfloatingiceEnum)) this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index])); - if (iomodel->Data(MaskVertexongroundediceEnum)) this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index])); - if (iomodel->numbernodetoelementconnectivity) this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,(IssmInt)iomodel->numbernodetoelementconnectivity[io_index])); - if (analysis_type==DiagnosticHorizAnalysisEnum) this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index])); + if (iomodel->Data(MeshVertexonbedEnum)) + this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index])); + if (iomodel->Data(MeshVertexonsurfaceEnum)) + this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index])); + if (iomodel->Data(MaskVertexonfloatingiceEnum)) + this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index])); + if (iomodel->Data(MaskVertexongroundediceEnum)) + this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index])); + if (iomodel->numbernodetoelementconnectivity) + this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,(IssmInt)iomodel->numbernodetoelementconnectivity[io_index])); + if (analysis_type==DiagnosticHorizAnalysisEnum) + this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index])); /*set single point constraints: */ @@ -72,9 +81,15 @@ /*Diagnostic Horiz*/ + #ifdef _HAVE_DIAGNOSTIC_ if (analysis_type==DiagnosticHorizAnalysisEnum){ + + /*Coordinate system provided, convert to coord_system matrix*/ + _assert_(iomodel->Data(DiagnosticReferentialEnum)); + XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6); + if (dim==3){ /*We have a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */ - if (!iomodel->Data(MeshVertexonbedEnum)) _error_("iomodel->nodeonbed is NULL"); - if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL"); + _assert_(iomodel->Data(MeshVertexonbedEnum)); + _assert_(iomodel->Data(FlowequationVertexEquationEnum)); if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){ for(k=1;k<=gsize;k++) this->FreezeDof(k); @@ -97,6 +112,6 @@ } } - - } + } + #endif /*Diagnostic Hutter*/ @@ -202,5 +217,5 @@ memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid); memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type); - memcpy(marshalled_dataset,&referential,9*sizeof(double));marshalled_dataset+=9*sizeof(double); + memcpy(marshalled_dataset,&coord_system,9*sizeof(double));marshalled_dataset+=9*sizeof(double); /*marshall objects: */ @@ -248,5 +263,5 @@ memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid); memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type); - memcpy(&referential,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double); + memcpy(&coord_system,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double); /*demarshall objects: */ @@ -507,4 +522,14 @@ } /*}}}*/ +#ifdef _HAVE_DIAGNOSTIC_ +/*FUNCTION Node::GetCoordinateSystem{{{1*/ +void Node::GetCoordinateSystem(double* coord_system_out){ + + /*Copy coord_system*/ + for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l]; + +} +/*}}}*/ +#endif /*FUNCTION Node::SetVertexDof {{{1*/ void Node::SetVertexDof(int in_dof){ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 10366) +++ /issm/trunk/src/c/objects/Node.h (revision 10367) @@ -30,5 +30,5 @@ Inputs* inputs; //properties of this node int analysis_type; - double referential[3][3]; + double coord_system[3][3]; /*Node constructors, destructors {{{1*/ @@ -65,38 +65,41 @@ /*}}}*/ /*Node numerical routines {{{1*/ - void Configure(DataSet* nodes,Vertices* vertices); - void CreateNodalConstraints(Vec ys); - void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); - int Sid(void); - int GetVertexDof(void); - int GetVertexId(void); - void SetVertexDof(int in_dof); - bool InAnalysis(int analysis_type); - int GetApproximation(); - int GetNumberOfDofs(int approximation_enum,int setenum); - int IsClone(); - void ApplyConstraint(int dof,double value); - void RelaxConstraint(int dof); - void DofInSSet(int dof); - void DofInFSet(int dof); - int GetDof(int dofindex,int setenum); - void CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s); - int GetConnectivity(); - void GetDofList(int* poutdoflist,int approximation_enum,int setenum); - void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); - int GetDofList1(void); - int GetSidList(void); + void Configure(DataSet* nodes,Vertices* vertices); + void CreateNodalConstraints(Vec ys); + void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); + int Sid(void); + int GetVertexDof(void); + int GetVertexId(void); +#ifdef _HAVE_DIAGNOSTIC_ + void GetCoordinateSystem(double* coord_system_out); +#endif + void SetVertexDof(int in_dof); + bool InAnalysis(int analysis_type); + int GetApproximation(); + int GetNumberOfDofs(int approximation_enum,int setenum); + int IsClone(); + void ApplyConstraint(int dof,double value); + void RelaxConstraint(int dof); + void DofInSSet(int dof); + void DofInFSet(int dof); + int GetDof(int dofindex,int setenum); + void CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s); + int GetConnectivity(); + void GetDofList(int* poutdoflist,int approximation_enum,int setenum); + void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); + int GetDofList1(void); + int GetSidList(void); double GetX(); double GetY(); double GetZ(); double GetSigma(); - int IsOnBed(); - int IsOnSurface(); - void FreezeDof(int dof); - int IsFloating(); - int IsGrounded(); - void UpdateSpcs(double* ys); - void VecMerge(Vec ug, double* vector_serial,int setnum); - void VecReduce(Vec vector, double* ug_serial,int setnum); + int IsOnBed(); + int IsOnSurface(); + void FreezeDof(int dof); + int IsFloating(); + int IsGrounded(); + void UpdateSpcs(double* ys); + void VecMerge(Vec ug, double* vector_serial,int setnum); + void VecReduce(Vec vector, double* ug_serial,int setnum); /*}}}*/ Index: /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp =================================================================== --- /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp (revision 10367) +++ /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp (revision 10367) @@ -0,0 +1,53 @@ +/*!\file: GetGlobalDofList.cpp + * \brief create transform matrix for different coordinate systems + */ +#include "./elements.h" + +void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension){ + + int i; + double *transform = NULL; + double *values = NULL; + double coord_system[3][3]; + + /*Some checks*/ + _assert_(numnodes && nodes); + if(dimension!=2 && dimension!=3) _error_("only 2d and 3d coordinate systems are supported"); + + /*Allocate and initialize transform matrix*/ + transform=(double*)xmalloc(numnodes*dimension*numnodes*dimension*sizeof(double)); + for(i=0;iGetCoordinateSystem(&coord_system[0][0]); + if(dimension==2){ + transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0]; + transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1]; + transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0]; + transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1]; + } + else if(dimension==3){ + transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0]; + transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1]; + transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+2] = coord_system[0][2]; + transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0]; + transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1]; + transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+2] = coord_system[1][2]; + transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+0] = coord_system[2][0]; + transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+1] = coord_system[2][1]; + transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+2] = coord_system[2][2]; + } + } + + /*Assign output pointer*/ + *ptransform=transform; +} Index: /issm/trunk/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk/src/c/shared/Elements/elements.h (revision 10366) +++ /issm/trunk/src/c/shared/Elements/elements.h (revision 10367) @@ -17,4 +17,5 @@ int* GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum); int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation_enum); +void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension); inline void printarray(double* array,int lines,int cols=1){ Index: /issm/trunk/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp =================================================================== --- /issm/trunk/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp (revision 10367) +++ /issm/trunk/src/c/shared/Numerics/XZvectorsToCoordinateSystem.cpp (revision 10367) @@ -0,0 +1,53 @@ +#include "../Matlab/matlabshared.h" +#include "../Alloc/alloc.h" +#include "../../include/include.h" +#include "../Exceptions/exceptions.h" +#include + +void XZvectorsToCoordinateSystem(double* T,double* xzvectors){ + + int i,j; + double x[3],y[3],z[3]; + double x_norm, y_norm, z_norm; + + for(i=0;i<6;i++){ + if(isnan(xzvectors[i])){ + /*At least one NaN found: default to Id*/ + T[0*3+0] = 1.0; T[0*3+1] = 0.0; T[0*3+2] = 0.0; + T[1*3+0] = 0.0; T[1*3+1] = 1.0; T[1*3+2] = 0.0; + T[2*3+0] = 0.0; T[2*3+1] = 0.0; T[2*3+2] = 1.0; + + return; + } + } + + /* get input {x} (vector in local x-z plane): */ + x[0] = xzvectors[0]; + x[1] = xzvectors[1]; + x[2] = xzvectors[2]; + + /* get input {z} (local tangent plane normal vector): */ + z[0] = xzvectors[3]; + z[1] = xzvectors[4]; + z[2] = xzvectors[5]; + + /* compute {y} = {z} x {x}: */ + y[0] = x[2]*z[1] - x[1]*z[2]; + y[1] = -x[2]*z[0] + x[0]*z[2]; + y[2] = x[1]*z[0] - x[0]*z[1]; + + /* normalise {x}, {y} and {z} to form unit vectors {i_hat}, {j_hat} and {k_hat}; + store in {x}, {y}, and {z}: */ + x_norm = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2]); + y_norm = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2]); + z_norm = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2]); + + x[0] = x[0]/x_norm; x[1] = x[1]/x_norm; x[2] = x[2]/x_norm; + y[0] = y[0]/y_norm; y[1] = y[1]/y_norm; y[2] = y[2]/y_norm; + z[0] = z[0]/z_norm; z[1] = z[1]/z_norm; z[2] = z[2]/z_norm; + + /* Tlg columns are just {i_hat}, {j_hat} and {k_hat}, respectively: */ + T[0*3+0] = x[0]; T[0*3+1] = y[0]; T[0*3+2] = z[0]; + T[1*3+0] = x[1]; T[1*3+1] = y[1]; T[1*3+2] = z[1]; + T[2*3+0] = x[2]; T[2*3+1] = y[2]; T[2*3+2] = z[2]; +} Index: /issm/trunk/src/c/shared/Numerics/numerics.h =================================================================== --- /issm/trunk/src/c/shared/Numerics/numerics.h (revision 10366) +++ /issm/trunk/src/c/shared/Numerics/numerics.h (revision 10367) @@ -28,4 +28,5 @@ double UnitConversion(double value, int direction_enum, int type_enum); void PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type); +void XZvectorsToCoordinateSystem(double* T,double* xzvectors); #endif //ifndef _NUMERICS_H_ Index: /issm/trunk/src/c/solutions/issm.cpp =================================================================== --- /issm/trunk/src/c/solutions/issm.cpp (revision 10366) +++ /issm/trunk/src/c/solutions/issm.cpp (revision 10367) @@ -107,5 +107,6 @@ #ifdef _HAVE_CONTROL_ #ifdef _HAVE_TAO_ - controltao_core(femmodel); + //controltao_core(femmodel); + control_core(femmodel); #else control_core(femmodel);