Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3354) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3355) @@ -22,4 +22,6 @@ using namespace std; +/*Constructors/Destructors*/ +/*FUNCTION DataSet::DataSet(){{{1*/ DataSet::DataSet(){ @@ -29,5 +31,6 @@ } - +/*}}}*/ +/*FUNCTION DataSet::DataSet(int dataset_enum){{{1*/ DataSet::DataSet(int dataset_enum){ enum_type=dataset_enum; @@ -38,5 +41,37 @@ } - +/*}}}*/ +/*FUNCTION DataSet::Copy{{{1*/ +DataSet* DataSet::Copy(void){ + + + DataSet* copy=NULL; + vector::iterator object; + Object* object_copy=NULL; + + copy=new DataSet(enum_type); + + copy->sorted=sorted; + copy->presorted=presorted; + if(sorted_ids){ + copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int)); + memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int)); + } + if(id_offsets){ + copy->id_offsets=(int*)xmalloc(objects.size()*sizeof(int)); + memcpy(copy->id_offsets,id_offsets,objects.size()*sizeof(int)); + } + + /*Now we need to deep copy the objects: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Call echo on object: */ + object_copy = (*object)->copy(); + copy->AddObject(object_copy); + } + return copy; +} +/*}}}*/ +/*FUNCTION DataSet::~DataSet{{{1*/ DataSet::~DataSet(){ clear(); @@ -44,65 +79,8 @@ xfree((void**)&id_offsets); } - -int DataSet::GetEnum(){ - return enum_type; -} - -void DataSet::Echo(){ - - - vector::iterator object; - - if(this==NULL)ISSMERROR(" trying to echo a NULL dataset"); - - _printf_("DataSet echo: %i objects\n",objects.size()); - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Call echo on object: */ - (*object)->Echo(); - - } - return; -} - -void DataSet::DeepEcho(){ - - - vector::iterator object; - - if(this==NULL)ISSMERROR(" trying to echo a NULL dataset"); - - _printf_("DataSet echo: %i objects\n",objects.size()); - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Call deep echo on object: */ - (*object)->DeepEcho(); - - } - return; -} - -int DataSet::MarshallSize(){ - - vector::iterator object; - int marshalled_dataset_size=0; - - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - marshalled_dataset_size+= (*object)->MarshallSize(); - } - - marshalled_dataset_size+=sizeof(int); //objects size - marshalled_dataset_size+=sizeof(int); //sorted size - if(sorted){ - marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids - marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets - } - - return marshalled_dataset_size; -} - +/*}}}*/ + +/*I/O*/ +/*FUNCTION DataSet::Marshall{{{1*/ char* DataSet::Marshall(){ @@ -150,33 +128,31 @@ return marshalled_dataset; } - -int DataSet::AddObject(Object* object){ - - objects.push_back(object); - - return 1; -} - -int DataSet::DeleteObject(Object* object){ - - vector::iterator iterator; - - if(object){ - iterator = find(objects.begin(), objects.end(),object); - delete *iterator; - objects.erase(iterator); - } - -} - -int DataSet::Size(void){ - - return objects.size(); -} - +/*}}}*/ +/*FUNCTION DataSet::MarshallSize{{{1*/ +int DataSet::MarshallSize(){ + + vector::iterator object; + int marshalled_dataset_size=0; + + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + marshalled_dataset_size+= (*object)->MarshallSize(); + } + + marshalled_dataset_size+=sizeof(int); //objects size + marshalled_dataset_size+=sizeof(int); //sorted size + if(sorted){ + marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids + marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets + } + + return marshalled_dataset_size; +} +/*}}}*/ +/*FUNCTION DataSetDemarshall{{{1*/ DataSet* DataSetDemarshall(char* marshalled_dataset){ int i; - + DataSet* dataset=NULL; int numobjects=0; @@ -203,10 +179,10 @@ } - #ifdef _ISSM_DEBUG_ +#ifdef _ISSM_DEBUG_ _printf_("Number of objects in dataset being demarshalled: %i\n",numobjects); - #endif +#endif for(i=0;i::iterator object; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + delete (*object); + } + objects.clear(); +} +/*}}}*/ +/*FUNCTION DataSet::DeleteObject{{{1*/ +int DataSet::DeleteObject(Object* object){ + + vector::iterator iterator; + + if(object){ + iterator = find(objects.begin(), objects.end(),object); + delete *iterator; + objects.erase(iterator); + } + +} +/*}}}*/ +/*FUNCTION DataSet::DeepEcho{{{1*/ +void DataSet::DeepEcho(){ + + + vector::iterator object; + + if(this==NULL)ISSMERROR(" trying to echo a NULL dataset"); + + _printf_("DataSet echo: %i objects\n",objects.size()); + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Call deep echo on object: */ + (*object)->DeepEcho(); + + } + return; +} +/*}}}*/ +/*FUNCTION DataSet::Echo{{{1*/ +void DataSet::Echo(){ + + + vector::iterator object; + + if(this==NULL)ISSMERROR(" trying to echo a NULL dataset"); + + _printf_("DataSet echo: %i objects\n",objects.size()); + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Call echo on object: */ + (*object)->Echo(); + + } + return; +} +/*}}}*/ +/*FUNCTION DataSet::FindParam(double* pscalar, char* name){{{1*/ int DataSet::FindParam(double* pscalar, char* name){ @@ -339,5 +388,6 @@ return found; } - +/*}}}*/ +/*FUNCTION DataSet::FindParam(int* pinteger,char* name){{{1*/ int DataSet::FindParam(int* pinteger,char* name){ @@ -369,5 +419,6 @@ return found; } - +/*}}}*/ +/*FUNCTION DataSet::FindParam(char** pstring,char* name){{{1*/ int DataSet::FindParam(char** pstring,char* name){ @@ -399,4 +450,6 @@ } +/*}}}*/ +/*FUNCTION DataSet::FindParam(char*** pstringarray,int* pM,char* name){{{1*/ int DataSet::FindParam(char*** pstringarray,int* pM,char* name){ @@ -429,4 +482,6 @@ } +/*}}}*/ +/*FUNCTION DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){{{1*/ int DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){ @@ -460,4 +515,6 @@ } +/*}}}*/ +/*FUNCTION DataSet::FindParam(Vec* pvec,char* name){{{1*/ int DataSet::FindParam(Vec* pvec,char* name){ @@ -489,4 +546,6 @@ } +/*}}}*/ +/*FUNCTION DataSet::FindParamMat* pmat,char* name){{{1*/ int DataSet::FindParam(Mat* pmat,char* name){ @@ -518,5 +577,32 @@ } - +/*}}}*/ +/*FUNCTION DataSet::FindParamObject{{{1*/ +Object* DataSet::FindParamObject(char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum()){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Return the object: */ + return (*object); + } + } + } + return NULL; +} +/*}}}*/ +/*FUNCTION DataSet::FindResult(Vec* presult,char* name){{{1*/ int DataSet::FindResult(Vec* presult,char* name){ @@ -547,5 +633,6 @@ return found; } - +/*}}}*/ +/*FUNCTION DataSet::FindResult(void* pvalue, char* name){{{1*/ int DataSet::FindResult(void* pvalue, char* name){ @@ -578,49 +665,66 @@ return found; } - - -Object* DataSet::FindParamObject(char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum()){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Return the object: */ - return (*object); - } - } - } - return NULL; -} - - +/*}}}*/ +/*FUNCTION DataSet::GetEnum(){{{1*/ +int DataSet::GetEnum(){ + return enum_type; +} +/*}}}*/ +/*FUNCTION DataSet::GetEnum(int offset){{{1*/ +int DataSet::GetEnum(int offset){ + + return objects[offset]->Enum(); + +} +/*}}}*/ +/*FUNCTION DataSet::GetObjectByOffset{{{1*/ +Object* DataSet::GetObjectByOffset(int offset){ + + return objects[offset]; + +} +/*}}}*/ +/*FUNCTION DataSet::GetObjectById{{{1*/ +Object* DataSet::GetObjectById(int* poffset,int eid){ + + int id_offset; + int offset; + int i; + + if(!sorted)ISSMERROR(" trying to binary search on a non-sorted dataset!"); + + /*Carry out a binary search on the sorted_ids: */ + if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){ + ISSMERROR(exprintf("%s%i"," could not find object with id ",eid)); + } + + /*Convert the id offset into sorted offset: */ + offset=id_offsets[id_offset]; + + /*Assign output pointers if requested:*/ + if (poffset)*poffset=offset; + + /*Return object at offset position in objects :*/ + return objects[offset]; +} +/*}}}*/ +/*FUNCTION DataSet::NodeRank{{{1*/ void DataSet::NodeRank(int* ranks){ /*Go through a dataset, and find a Node* object. For this object, * ask it to report its cpu rank: */ - + int node_rank; int id; vector::iterator object; - + for ( object=objects.begin() ; object < objects.end(); object++ ){ /*Call echo on object: */ if((*object)->Enum()==NodeEnum()){ - + /*Ok, this object is a node, ask which rank it has: */ node_rank=(*object)->MyRank(); - + /*Which id does it have: */ id=(*object)->GetId(); @@ -632,6 +736,213 @@ return; } - - +/*}}}*/ +/*FUNCTION DataSet::Presort{{{1*/ +void DataSet::Presort(){ + + /*vector of objects is already sorted, just allocate the sorted ids and their + * offsets:*/ + int i; + + if(objects.size()){ + sorted_ids=(int*)xmalloc(objects.size()*sizeof(int)); + id_offsets=(int*)xmalloc(objects.size()*sizeof(int)); + for(i=0;iGetId(); + } + } + + /*set sorted flag: */ + sorted=1; +} +/*}}}*/ +/*FUNCTION DataSet::SetSorting{{{1*/ +void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){ + + sorted=1; + sorted_ids=in_sorted_ids; + id_offsets=in_id_offsets; +} +/*}}}*/ +/*FUNCTION DataSet::Size{{{1*/ +int DataSet::Size(void){ + + return objects.size(); +} +/*}}}*/ +/*FUNCTION DataSet::Sort{{{1*/ +void DataSet::Sort(){ + + /*Only sort if we are not already sorted: */ + if(!sorted){ + ISSMERROR(" not implemented yet!"); + } +} +/*}}}*/ + +/*Objects methods*/ +/*FUNCTION DataSet::ComputePressure{{{1*/ +void DataSet::ComputePressure(Vec p_g){ + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->ComputePressure(p_g); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::Configure{{{1*/ +void DataSet::Configure(DataSet* elements,DataSet* loads, DataSet* nodes, DataSet* materials,DataSet* parameters){ + + vector::iterator object; + Element* element=NULL; + Load* load=NULL; + Node* node=NULL; + Numpar* numpar=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->Configure(loads,nodes,materials,parameters); + } + if(EnumIsLoad((*object)->Enum())){ + + load=(Load*)(*object); + + load->Configure(elements,nodes,materials); + } + + if((*object)->Enum()==NodeEnum()){ + node=(Node*)(*object); + node->Configure(nodes); + } + + if((*object)->Enum()==NumparEnum()){ + numpar=(Numpar*)(*object); + numpar->Configure(parameters); + } + + } + +} +/*}}}*/ +/*FUNCTION DataSet::CostFunction{{{1*/ +void DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ + + double J=0;; + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + J+=element->CostFunction(inputs,analysis_type,sub_analysis_type); + + } + } + + /*Assign output pointers:*/ + *pJ=J; + +} +/*}}}*/ +/*FUNCTION DataSet::CreateKMatrix{{{1*/ +void DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + + vector::iterator object; + Element* element=NULL; + Load* load=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); + } + if(EnumIsLoad((*object)->Enum())){ + + load=(Load*)(*object); + load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::CreatePartitioningVector{{{1*/ +void DataSet::CreatePartitioningVector(Vec* ppartition,int numberofnodes,int numdofspernode){ + + /*output: */ + Vec partition=NULL; + vector::iterator object; + Node* node=NULL; + + /*Create partition vector: */ + partition=NewVec(numberofnodes); + + /*Go through all nodes, and ask each node to plug its 1D doflist in + * partition. The location where each node plugs its doflist into + * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs). + */ + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Check this is a node: */ + if((*object)->Enum()==NodeEnum()){ + + node=(Node*)(*object); + + /*Ok, this object is a node, ask it to plug values into partition: */ + node->CreatePartition(partition); + + } + } + + /*Assemble the petsc vector: */ + VecAssemblyBegin(partition); + VecAssemblyEnd(partition); + + /*Assign output pointers: */ + *ppartition=partition; + + return; +} +/*}}}*/ +/*FUNCTION DataSet::CreatePVector{{{1*/ +void DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ + + vector::iterator object; + Element* element=NULL; + Load* load=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); + } + if(EnumIsLoad((*object)->Enum())){ + + load=(Load*)(*object); + load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::DistributeDofs{{{1*/ void DataSet::DistributeDofs(int numberofnodes,int numdofspernode){ @@ -650,5 +961,5 @@ extern int num_procs; extern int my_rank; - + dofcount=0; for ( object=objects.begin() ; object < objects.end(); object++ ){ @@ -658,8 +969,8 @@ node=(Node*)(*object); - + /*Ok, this object is a node, ask it to distribute dofs, and update the dofcount: */ node->DistributeDofs(&dofcount,&dofcount1); - + } } @@ -668,5 +979,5 @@ *0. This means the dofs between all the cpus are not synchronized! We need to synchronize all *dof on all cpus, and use those to update the dofs of every node: */ - + alldofcount=(int*)xmalloc(num_procs*sizeof(int)); MPI_Gather(&dofcount,1,MPI_INT,alldofcount,1,MPI_INT,0,MPI_COMM_WORLD); @@ -700,8 +1011,8 @@ node=(Node*)(*object); - + /*Ok, this object is a node, ask it to update his dofs: */ node->UpdateDofs(dofcount,dofcount1); - + } } @@ -720,8 +1031,8 @@ node=(Node*)(*object); - + /*Ok, let this object show its border dofs, if is is a border dof: */ node->ShowBorderDofs(borderdofs,borderdofs1); - + } } @@ -736,8 +1047,8 @@ node=(Node*)(*object); - + /*Ok, this object is a node, ask it to update his dofs: */ node->UpdateBorderDofs(allborderdofs,allborderdofs1); - + } } @@ -755,49 +1066,72 @@ } - - -void DataSet::CreatePartitioningVector(Vec* ppartition,int numberofnodes,int numdofspernode){ - - /*output: */ - Vec partition=NULL; +/*}}}*/ +/*FUNCTION DataSet::Du{{{1*/ +void DataSet::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){ + + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->Du(du_g,inputs,analysis_type,sub_analysis_type); + } + } + + +} +/*}}}*/ +/*FUNCTION DataSet::FieldDepthAverageAtBase{{{1*/ +void DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){ + vector::iterator object; Node* node=NULL; - /*Create partition vector: */ - partition=NewVec(numberofnodes); - - /*Go through all nodes, and ask each node to plug its 1D doflist in - * partition. The location where each node plugs its doflist into - * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs). - */ - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Check this is a node: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + if((*object)->Enum()==NodeEnum()){ - + node=(Node*)(*object); - - /*Ok, this object is a node, ask it to plug values into partition: */ - node->CreatePartition(partition); - - } - } - - /*Assemble the petsc vector: */ - VecAssemblyBegin(partition); - VecAssemblyEnd(partition); - - /*Assign output pointers: */ - *ppartition=partition; - - return; -} - + node->FieldDepthAverageAtBase(field,field_serial,fieldname); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::FieldExtrude{{{1*/ +void DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){ + + vector::iterator object; + Penta* penta=NULL; + Node* node=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum()){ + + penta=(Penta*)(*object); + penta->FieldExtrude(field,field_serial,field_name,collapse); + + } + if((*object)->Enum()==NodeEnum()){ + + node=(Node*)(*object); + node->FieldExtrude(field,field_serial,field_name); + + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::FlagClones{{{1*/ void DataSet::FlagClones(int numberofnodes){ int i; extern int num_procs; - + int* ranks=NULL; int* minranks=NULL; @@ -805,9 +1139,9 @@ vector::iterator object; Node* node=NULL; - + /*Allocate ranks: */ ranks=(int*)xmalloc(numberofnodes*sizeof(int)); minranks=(int*)xmalloc(numberofnodes*sizeof(int)); - + for(i=0;iEnum()==NodeEnum()){ - + node=(Node*)(*object); - + /*Ok, this object is a node, ask it to plug values into partition: */ node->SetClone(minranks); - + } } @@ -852,6 +1186,152 @@ } +/*}}}*/ +/*FUNCTION DataSet::FlagNodeSets{{{1*/ +void DataSet::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s){ + + vector::iterator object; + Node* node=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Check this is a single point constraint (spc): */ + if((*object)->Enum()==NodeEnum()){ + + node=(Node*)(*object); + + /*Plug set values intp our 4 set vectors: */ + node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s); + + } + } + + /*Assemble: */ + VecAssemblyBegin(pv_g); + VecAssemblyEnd(pv_g); + + VecAssemblyBegin(pv_m); + VecAssemblyEnd(pv_m); + + VecAssemblyBegin(pv_n); + VecAssemblyEnd(pv_n); + + VecAssemblyBegin(pv_f); + VecAssemblyEnd(pv_f); + + VecAssemblyBegin(pv_s); + VecAssemblyEnd(pv_s); + +} +/*}}}*/ +/*FUNCTION DataSet::Gradj{{{1*/ +void DataSet::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ + + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->Gradj(grad_g,inputs,analysis_type,sub_analysis_type,control_type); + } + } + + +} +/*}}}*/ +/*FUNCTION DataSet::MeltingIsPresent{{{1*/ +int DataSet::MeltingIsPresent(){ + + int found=0; + int mpi_found=0; + + vector::iterator object; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PengridEnum()){ + found=1; + break; + } + } + + #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::MeltingConstraints{{{1*/ +void DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ + + /* generic object pointer: */ + vector::iterator object; + Pengrid* pengrid=NULL; + + int unstable=0; + int num_unstable_constraints=0; + int converged=0; + int sum_num_unstable_constraints=0; + + num_unstable_constraints=0; + + /*Enforce constraints: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ - + if((*object)->Enum()==PengridEnum()){ + + pengrid=(Pengrid*)(*object); + + pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Have we converged? : */ + if (num_unstable_constraints==0) converged=1; + else converged=0; + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} +/*}}}*/ +/*FUNCTION DataSet::Misfit{{{1*/ +void DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ + + double J=0;; + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + J+=element->Misfit(inputs,analysis_type,sub_analysis_type); + + } + } + + /*Assign output pointers:*/ + *pJ=J; + +} +/*}}}*/ +/*FUNCTION DataSet::NumberOfDofs{{{1*/ int DataSet::NumberOfDofs(){ @@ -867,12 +1347,12 @@ /*Check this is a node: */ if((*object)->Enum()==NodeEnum()){ - + node=(Node*)(*object); - + /*Ok, this object is a node, ask it to plug values into partition: */ if (!node->IsClone()){ numdofs+=node->GetNumberOfDofs(); - + } } @@ -884,45 +1364,6 @@ return allnumdofs; } - -void DataSet::SetupSpcs(DataSet* nodes,Vec yg){ - - 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); - - /*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::NumberOfRgbs{{{1*/ int DataSet::NumberOfRgbs(){ @@ -943,5 +1384,87 @@ return count; } - +/*}}}*/ +/*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::PenaltyCreateKMatrix{{{1*/ +void DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ + + vector::iterator object; + Load* load=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsLoad((*object)->Enum())){ + + load=(Load*)(*object); + load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::PenaltyCreatePVector{{{1*/ +void DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ + + vector::iterator object; + Load* load=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsLoad((*object)->Enum())){ + + load=(Load*)(*object); + load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type); + } + } + +} +/*}}}*/ +/*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::SetupMpcs{{{1*/ void DataSet::SetupMpcs(Mat Rmg,DataSet* nodes){ @@ -997,5 +1520,5 @@ node1->DofInMSet(dof-1); } - + /*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the *same: */ @@ -1010,195 +1533,68 @@ } } - -void DataSet::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s){ - - vector::iterator object; +/*}}}*/ +/*FUNCTION DataSet::SetupSpcs{{{1*/ +void DataSet::SetupSpcs(DataSet* nodes,Vec yg){ + + 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()==NodeEnum()){ - - node=(Node*)(*object); - - /*Plug set values intp our 4 set vectors: */ - node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s); - - } - } - - /*Assemble: */ - VecAssemblyBegin(pv_g); - VecAssemblyEnd(pv_g); - - VecAssemblyBegin(pv_m); - VecAssemblyEnd(pv_m); - - VecAssemblyBegin(pv_n); - VecAssemblyEnd(pv_n); - - VecAssemblyBegin(pv_f); - VecAssemblyEnd(pv_f); - - VecAssemblyBegin(pv_s); - VecAssemblyEnd(pv_s); - -} -void DataSet::clear(){ - - vector::iterator object; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - delete (*object); - } - objects.clear(); -} - - -void DataSet::Configure(DataSet* elements,DataSet* loads, DataSet* nodes, DataSet* materials,DataSet* parameters){ + if((*object)->Enum()==SpcEnum()){ + + spc=(Spc*)(*object); + + /*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::SurfaceArea{{{1*/ +void DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){ + + double S=0;; vector::iterator object; Element* element=NULL; - Load* load=NULL; - Node* node=NULL; - Numpar* numpar=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + if(EnumIsElement((*object)->Enum())){ element=(Element*)(*object); - element->Configure(loads,nodes,materials,parameters); - } - if(EnumIsLoad((*object)->Enum())){ - - load=(Load*)(*object); - - load->Configure(elements,nodes,materials); - } - - if((*object)->Enum()==NodeEnum()){ - node=(Node*)(*object); - node->Configure(nodes); - } - - if((*object)->Enum()==NumparEnum()){ - numpar=(Numpar*)(*object); - numpar->Configure(parameters); - } - - } - -} - -void DataSet::Presort(){ - - /*vector of objects is already sorted, just allocate the sorted ids and their - * offsets:*/ - int i; - - if(objects.size()){ - sorted_ids=(int*)xmalloc(objects.size()*sizeof(int)); - id_offsets=(int*)xmalloc(objects.size()*sizeof(int)); - for(i=0;iGetId(); - } - } - - /*set sorted flag: */ - sorted=1; -} - -void DataSet::Sort(){ - - /*Only sort if we are not already sorted: */ - if(!sorted){ - ISSMERROR(" not implemented yet!"); - } -} - - -Object* DataSet::GetObjectByOffset(int offset){ - - return objects[offset]; - -} - -Object* DataSet::GetObjectById(int* poffset,int eid){ - - int id_offset; - int offset; - int i; - - if(!sorted)ISSMERROR(" trying to binary search on a non-sorted dataset!"); - - /*Carry out a binary search on the sorted_ids: */ - if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){ - ISSMERROR(exprintf("%s%i"," could not find object with id ",eid)); - } - - /*Convert the id offset into sorted offset: */ - offset=id_offsets[id_offset]; - - /*Assign output pointers if requested:*/ - if (poffset)*poffset=offset; - - /*Return object at offset position in objects :*/ - return objects[offset]; -} - -void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){ - - sorted=1; - sorted_ids=in_sorted_ids; - id_offsets=in_id_offsets; -} - -void DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - - vector::iterator object; - Element* element=NULL; - Load* load=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); - } - if(EnumIsLoad((*object)->Enum())){ - - load=(Load*)(*object); - load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type); - } - } - -} - -void DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ - - vector::iterator object; - Element* element=NULL; - Load* load=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); - } - if(EnumIsLoad((*object)->Enum())){ - - load=(Load*)(*object); - load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type); - } - } - -} - + S+=element->SurfaceArea(inputs,analysis_type,sub_analysis_type); + + } + } + + /*Assign output pointers:*/ + *pS=S; + +} +/*}}}*/ +/*FUNCTION DataSet::UpdateFromInputs{{{1*/ void DataSet::UpdateFromInputs(void* inputs){ @@ -1233,319 +1629,6 @@ } - -void DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ - - vector::iterator object; - Load* load=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsLoad((*object)->Enum())){ - - load=(Load*)(*object); - load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type); - } - } - -} - -void DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ - - vector::iterator object; - Load* load=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsLoad((*object)->Enum())){ - - load=(Load*)(*object); - load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type); - } - } - -} - -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; -} - -int DataSet::MeltingIsPresent(){ - - int found=0; - int mpi_found=0; - - vector::iterator object; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PengridEnum()){ - found=1; - break; - } - } - - #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; -} - -void DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ - - /* generic object pointer: */ - vector::iterator object; - Pengrid* pengrid=NULL; - - int unstable=0; - int num_unstable_constraints=0; - int converged=0; - int sum_num_unstable_constraints=0; - - num_unstable_constraints=0; - - /*Enforce constraints: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PengridEnum()){ - - pengrid=(Pengrid*)(*object); - - pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type); - - num_unstable_constraints+=unstable; - } - } - - #ifdef _PARALLEL_ - MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); - MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); - num_unstable_constraints=sum_num_unstable_constraints; - #endif - - /*Have we converged? : */ - if (num_unstable_constraints==0) converged=1; - else converged=0; - - /*Assign output pointers: */ - *pconverged=converged; - *pnum_unstable_constraints=num_unstable_constraints; -} - -DataSet* DataSet::Copy(void){ - - - DataSet* copy=NULL; - vector::iterator object; - Object* object_copy=NULL; - - copy=new DataSet(enum_type); - - copy->sorted=sorted; - copy->presorted=presorted; - if(sorted_ids){ - copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int)); - memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int)); - } - if(id_offsets){ - copy->id_offsets=(int*)xmalloc(objects.size()*sizeof(int)); - memcpy(copy->id_offsets,id_offsets,objects.size()*sizeof(int)); - } - - /*Now we need to deep copy the objects: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Call echo on object: */ - object_copy = (*object)->copy(); - copy->AddObject(object_copy); - } - return copy; -} - - -void DataSet::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){ - - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - element->Du(du_g,inputs,analysis_type,sub_analysis_type); - } - } - - -} - -void DataSet::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ - - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - element->Gradj(grad_g,inputs,analysis_type,sub_analysis_type,control_type); - } - } - - -} - -void DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ - - double J=0;; - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - J+=element->Misfit(inputs,analysis_type,sub_analysis_type); - - } - } - - /*Assign output pointers:*/ - *pJ=J; - -} - -void DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){ - - double J=0;; - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - J+=element->CostFunction(inputs,analysis_type,sub_analysis_type); - - } - } - - /*Assign output pointers:*/ - *pJ=J; - -} - -void DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){ - - double S=0;; - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - S+=element->SurfaceArea(inputs,analysis_type,sub_analysis_type); - - } - } - - /*Assign output pointers:*/ - *pS=S; - -} - -void DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){ - - vector::iterator object; - Penta* penta=NULL; - Node* node=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PentaEnum()){ - - penta=(Penta*)(*object); - penta->FieldExtrude(field,field_serial,field_name,collapse); - - } - if((*object)->Enum()==NodeEnum()){ - - node=(Node*)(*object); - node->FieldExtrude(field,field_serial,field_name); - - } - } - -} - -void DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){ - - vector::iterator object; - Node* node=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==NodeEnum()){ - - node=(Node*)(*object); - node->FieldDepthAverageAtBase(field,field_serial,fieldname); - } - } - -} - -void DataSet::ComputePressure(Vec p_g){ - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if(EnumIsElement((*object)->Enum())){ - - element=(Element*)(*object); - element->ComputePressure(p_g); - } - } - -} - - +/*}}}*/ +/*FUNCTION DataSet::UpdateNodePositions{{{1*/ void DataSet::UpdateNodePositions(double* thickness,double* bed){ @@ -1562,27 +1645,3 @@ } } - - -int DataSet::GetEnum(int offset){ - - return objects[offset]->Enum(); - -} - -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); - } - } - - - -} +/*}}}*/