Index: /issm/trunk/src/c/DataSet/Constraints.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Constraints.cpp (revision 4219) +++ /issm/trunk/src/c/DataSet/Constraints.cpp (revision 4220) @@ -64,3 +64,139 @@ } /*}}}*/ - +/*FUNCTION Constraints::NumberOfConstraints{{{1*/ +int Constraints::NumberOfConstraints(void){ + + int localconstraints; + int numberofconstraints; + + /*Get number of local constraints*/ + localconstraints=this->Size(); + + /*figure out total number of constraints combining all the cpus (no clones here)*/ + #ifdef _PARALLEL_ + MPI_Reduce(&localconstraints,&numberofconstraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&numberofconstraints,1,MPI_INT,0,MPI_COMM_WORLD); + #else + numberofconstraints=localconstraints; + #endif + + return numberofconstraints; +} +/*}}}*/ +/*FUNCTION Constraints::SetupSpcs{{{1*/ +void Constraints::SetupSpcs(Nodes* nodes,Vec yg,int analysis_type){ + + int i; + + Node* node=NULL; + int nodeid; + int dof; + double value; + + for(i=0;iSize();i++){ + + Object* object=(Object*)this->GetObjectByOffset(i); + + /*Check this is a single point constraint (spc): */ + if(object->Enum()==SpcEnum){ + + Spc* spc=(Spc*)object; + + if(spc->InAnalysis(analysis_type)){ + + /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */ + nodeid=spc->GetNodeId(); + dof=spc->GetDof(); + value=spc->GetValue(); + + /*Now, chase through nodes and find the corect node: */ + node=(Node*)nodes->GetObjectById(NULL,nodeid); + + /*Apply constraint: */ + if(node){ //in case the spc is dealing with a node on another cpu + node->ApplyConstraint(yg,dof,value); + } + } + + } + } + + /*Assemble yg: */ + VecAssemblyBegin(yg); + VecAssemblyEnd(yg); +} +/*}}}*/ +/*FUNCTION Constraints::SetupMpcs{{{1*/ +void Constraints::SetupMpcs(Mat Rmg,Nodes* nodes,int analysis_type){ + + int i; + + int nodeid1; + int nodeid2; + int dof; + + int dof1; + int dof2; + + + Node* node1=NULL; + Node* node2=NULL; + + int count=-1; + + for(i=0;iSize();i++){ + + Object* object=(Object*)this->GetObjectByOffset(i); + + /*Check this is a mutiple point constraint (spc): */ + if(object->Enum()==RgbEnum){ + + Rgb* rgb=(Rgb*)object; + + if(rgb->InAnalysis(analysis_type)){ + + /*we found an rgb, increment counter, so that row index for Rmg is up to date: */ + count++; + + + nodeid1=rgb->GetNodeId1(); + nodeid2=rgb->GetNodeId2(); + dof=rgb->GetDof(); + + /*For this rgb, find the nodes that go with it: */ + node1=(Node*)nodes->GetObjectById(NULL,nodeid1); + node2=(Node*)nodes->GetObjectById(NULL,nodeid2); + + if ((node1 && !node2) || (!node1 && node2)){ + /*we are missing one node, not good!*/ + ISSMERROR("%s%p%s%p"," in Rgb, missing one node. node1: ",node1," node2: ",node2); + } + + if(!node1 && !node2){ + /*That's ok, this Rgb can't find those nodes, so leave them alone. They are probably not on this + * cpu!*/ + } + else{ + /*Ok, this cpu owns both nodes. Put dof for node1 into m set, unless it is already there, + * in which case node2 gets into the m set: */ + if(node1->DofIsInMSet(dof-1)){ + node2->DofInMSet(dof-1); + } + else{ + node1->DofInMSet(dof-1); + } + + /*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the + *same: */ + dof1=node1->GetDof(dof-1); //matlab indexing + dof2=node2->GetDof(dof-1); //matlab indexing + + MatSetValue(Rmg,count,dof1,1.0,INSERT_VALUES); + MatSetValue(Rmg,count,dof2,-1.0,INSERT_VALUES); + + } + } + } + } +} +/*}}}*/ Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4219) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4220) @@ -562,6 +562,4 @@ } /*}}}*/ - -/*Objects methods*/ /*FUNCTION DataSet::Configure{{{1*/ void DataSet::Configure(Elements* elements,Loads* loads, DataSet* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){ @@ -592,233 +590,2 @@ } /*}}}*/ -/*FUNCTION DataSet::NumberOfVertices{{{1*/ -int DataSet::NumberOfVertices(void){ - - vector::iterator object; - int max_sid=0; - int sid; - int vertex_max_sid; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==VertexEnum){ - - Vertex* vertex=(Vertex*)(*object); - sid=vertex->Sid(); - if (sid>max_sid)max_sid=sid; - } - } - - #ifdef _PARALLEL_ - MPI_Reduce (&max_sid,&vertex_max_sid,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD ); - MPI_Bcast(&vertex_max_sid,1,MPI_INT,0,MPI_COMM_WORLD); - max_sid=vertex_max_sid; - #endif - - /*sid starts at 0*/ - max_sid++; - - /*return:*/ - return max_sid; -} -/*}}}*/ -/*FUNCTION DataSet::NumberOfLoads{{{1*/ -int DataSet::NumberOfLoads(void){ - - int localloads; - int numberofloads; - - /*Get number of local loads*/ - localloads=this->Size(); - - /*figure out total number of loads combining all the cpus (no clones here)*/ -#ifdef _PARALLEL_ - MPI_Reduce(&localloads,&numberofloads,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&numberofloads,1,MPI_INT,0,MPI_COMM_WORLD); -#else - numberofloads=localloads; -#endif - - return numberofloads; -} -/*}}}*/ -/*FUNCTION DataSet::NumberOfConstraints{{{1*/ -int DataSet::NumberOfConstraints(void){ - - int localconstraints; - int numberofconstraints; - - /*Get number of local constraints*/ - localconstraints=this->Size(); - - /*figure out total number of constraints combining all the cpus (no clones here)*/ -#ifdef _PARALLEL_ - MPI_Reduce(&localconstraints,&numberofconstraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&numberofconstraints,1,MPI_INT,0,MPI_COMM_WORLD); -#else - numberofconstraints=localconstraints; -#endif - - return numberofconstraints; -} -/*}}}*/ -/*FUNCTION DataSet::OutputRifts{{{1*/ -void DataSet::OutputRifts(Vec riftproperties){ - - vector::iterator object; - Riftfront* riftfront=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==RiftfrontEnum){ - - riftfront=(Riftfront*)(*object); - riftfront->OutputProperties(riftproperties); - } - } - - - -} -/*}}}*/ -/*FUNCTION DataSet::RiftIsPresent{{{1*/ -int DataSet::RiftIsPresent(){ - - int i; - - vector::iterator object; - Penpair* penpair=NULL; - int found=0; - int mpi_found=0; - - /*go though loads, and figure out if one of the loads is a PenPair with numdof=2: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PenpairEnum){ - - penpair=(Penpair*)(*object); - } - } - -#ifdef _PARALLEL_ - MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); - found=mpi_found; -#endif - - return found; -} -/*}}}*/ -/*FUNCTION DataSet::SetupSpcs{{{1*/ -void DataSet::SetupSpcs(DataSet* nodes,Vec yg,int analysis_type){ - - vector::iterator object; - Spc* spc=NULL; - Node* node=NULL; - - int nodeid; - int dof; - double value; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Check this is a single point constraint (spc): */ - if((*object)->Enum()==SpcEnum){ - - spc=(Spc*)(*object); - - if(spc->InAnalysis(analysis_type)){ - - /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */ - nodeid=spc->GetNodeId(); - dof=spc->GetDof(); - value=spc->GetValue(); - - /*Now, chase through nodes and find the corect node: */ - node=(Node*)nodes->GetObjectById(NULL,nodeid); - - /*Apply constraint: */ - if(node){ //in case the spc is dealing with a node on another cpu - node->ApplyConstraint(yg,dof,value); - } - } - - } - } - - /*Assemble yg: */ - VecAssemblyBegin(yg); - VecAssemblyEnd(yg); -} -/*}}}*/ -/*FUNCTION DataSet::SetupMpcs{{{1*/ -void DataSet::SetupMpcs(Mat Rmg,DataSet* nodes,int analysis_type){ - - vector::iterator object; - Rgb* rgb=NULL; - int nodeid1; - int nodeid2; - int dof; - - int dof1; - int dof2; - - - Node* node1=NULL; - Node* node2=NULL; - - int count=-1; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Check this is a single point constraint (spc): */ - if((*object)->Enum()==RgbEnum){ //we assume uniqueness of all Rgbs, no error checking here. - - rgb=(Rgb*)(*object); - if(rgb->InAnalysis(analysis_type)){ - - /*we found an rgb, increment counter, so that row index for Rmg is up to date: */ - count++; - - - nodeid1=rgb->GetNodeId1(); - nodeid2=rgb->GetNodeId2(); - dof=rgb->GetDof(); - - /*For this rgb, find the nodes that go with it: */ - node1=(Node*)nodes->GetObjectById(NULL,nodeid1); - node2=(Node*)nodes->GetObjectById(NULL,nodeid2); - - if ((node1 && !node2) || (!node1 && node2)){ - /*we are missing one node, not good!*/ - ISSMERROR("%s%p%s%p"," in Rgb, missing one node. node1: ",node1," node2: ",node2); - } - - if(!node1 && !node2){ - /*That's ok, this Rgb can't find those nodes, so leave them alone. They are probably not on this - * cpu!*/ - } - else{ - /*Ok, this cpu owns both nodes. Put dof for node1 into m set, unless it is already there, - * in which case node2 gets into the m set: */ - if(node1->DofIsInMSet(dof-1)){ - node2->DofInMSet(dof-1); - } - else{ - node1->DofInMSet(dof-1); - } - - /*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the - *same: */ - dof1=node1->GetDof(dof-1); //matlab indexing - dof2=node2->GetDof(dof-1); //matlab indexing - - MatSetValue(Rmg,count,dof1,1.0,INSERT_VALUES); - MatSetValue(Rmg,count,dof2,-1.0,INSERT_VALUES); - - } - } - } - } -} -/*}}}*/ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 4219) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 4220) @@ -59,9 +59,4 @@ Object* FindParamObject(char* name); void Ranks(int* ranks); - int NumberOfVertices(void); - int NumberOfLoads(void); - int NumberOfConstraints(void); - void SetupSpcs(DataSet* nodes,Vec yg,int analysis_type); - void SetupMpcs(Mat Rmg,DataSet* nodes,int analysis_type); void clear(); void Configure(Elements* elements,Loads* loads, DataSet* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); @@ -71,8 +66,6 @@ void SetSorting(int* in_sorted_ids,int* in_id_offsets); void Sort(); - int RiftIsPresent(); DataSet* Copy(void); int DeleteObject(Object* object); - void OutputRifts(Vec riftproperties); /*}}}*/ @@ -166,4 +159,5 @@ void DistributeDofs(int numberofnodes,int numdofspernode); void FlagClones(int numberofnodes); + int NumberOfVertices(void); /*}}}*/ @@ -185,4 +179,6 @@ int MeltingIsPresent(); void MeltingConstraints(int* pconverged, int* pnum_unstable_constraints); + int NumberOfLoads(void); + void OutputRifts(Vec riftproperties); /*}}}*/ @@ -218,4 +214,7 @@ /*numerics: {{{1*/ int NumberOfLocalRgbs(int analysis_type); + int NumberOfConstraints(void); + void SetupSpcs(Nodes* nodes,Vec yg,int analysis_type); + void SetupMpcs(Mat Rmg,Nodes* nodes,int analysis_type); /*}}}*/ Index: /issm/trunk/src/c/DataSet/Loads.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Loads.cpp (revision 4219) +++ /issm/trunk/src/c/DataSet/Loads.cpp (revision 4220) @@ -105,2 +105,39 @@ } /*}}}*/ +/*FUNCTION Loads::NumberOfLoads{{{1*/ +int Loads::NumberOfLoads(void){ + + int localloads; + int numberofloads; + + /*Get number of local loads*/ + localloads=this->Size(); + + /*figure out total number of loads combining all the cpus (no clones here)*/ + #ifdef _PARALLEL_ + MPI_Reduce(&localloads,&numberofloads,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&numberofloads,1,MPI_INT,0,MPI_COMM_WORLD); + #else + numberofloads=localloads; + #endif + + return numberofloads; +} +/*}}}*/ +/*FUNCTION Loads::OutputRifts{{{1*/ +void Loads::OutputRifts(Vec riftproperties){ + + int i; + + for(i=0;iSize();i++){ + + Object* object=(Object*)this->GetObjectByOffset(i); + + if(object->Enum()==RiftfrontEnum){ + + Riftfront* riftfront=(Riftfront*)object; + riftfront->OutputProperties(riftproperties); + } + } +} +/*}}}*/ Index: /issm/trunk/src/c/DataSet/Vertices.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Vertices.cpp (revision 4219) +++ /issm/trunk/src/c/DataSet/Vertices.cpp (revision 4220) @@ -179,2 +179,31 @@ } /*}}}*/ +/*FUNCTION Vertices::NumberOfVertices{{{1*/ +int Vertices::NumberOfVertices(void){ + + int i; + + int max_sid=0; + int sid; + int vertex_max_sid; + + for(i=0;iSize();i++){ + + Vertex* vertex=(Vertex*)this->GetObjectByOffset(i); + sid=vertex->Sid(); + if (sid>max_sid)max_sid=sid; + } + + #ifdef _PARALLEL_ + MPI_Reduce (&max_sid,&vertex_max_sid,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&vertex_max_sid,1,MPI_INT,0,MPI_COMM_WORLD); + max_sid=vertex_max_sid; + #endif + + /*sid starts at 0*/ + max_sid++; + + /*return:*/ + return max_sid; +} +/*}}}*/