Changeset 3355


Ignore:
Timestamp:
03/31/10 07:41:18 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added folds (really needed here)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r3332 r3355  
    2222using namespace std;
    2323
     24/*Constructors/Destructors*/
     25/*FUNCTION DataSet::DataSet(){{{1*/
    2426DataSet::DataSet(){
    2527       
    … …  
    2931
    3032}
    31 
     33/*}}}*/
     34/*FUNCTION DataSet::DataSet(int dataset_enum){{{1*/
    3235DataSet::DataSet(int dataset_enum){
    3336        enum_type=dataset_enum;
    … …  
    3841
    3942}
    40 
     43/*}}}*/
     44/*FUNCTION DataSet::Copy{{{1*/
     45DataSet*   DataSet::Copy(void){
     46
     47
     48        DataSet* copy=NULL;
     49        vector<Object*>::iterator object;
     50        Object* object_copy=NULL;
     51
     52        copy=new DataSet(enum_type);
     53
     54        copy->sorted=sorted;
     55        copy->presorted=presorted;
     56        if(sorted_ids){
     57                copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
     58                memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int));
     59        }
     60        if(id_offsets){
     61                copy->id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
     62                memcpy(copy->id_offsets,id_offsets,objects.size()*sizeof(int));
     63        }
     64
     65        /*Now we need to deep copy the objects: */
     66        for ( object=objects.begin() ; object < objects.end(); object++ ){
     67
     68                /*Call echo on object: */
     69                object_copy = (*object)->copy();
     70                copy->AddObject(object_copy);
     71        }
     72        return copy;
     73}
     74/*}}}*/
     75/*FUNCTION DataSet::~DataSet{{{1*/
    4176DataSet::~DataSet(){
    4277        clear();
    … …  
    4479        xfree((void**)&id_offsets);
    4580}
    46 
    47 int  DataSet::GetEnum(){
    48         return enum_type;
    49 }
    50        
    51 void DataSet::Echo(){
    52 
    53        
    54         vector<Object*>::iterator object;
    55 
    56         if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
    57 
    58         _printf_("DataSet echo: %i objects\n",objects.size());
    59 
    60         for ( object=objects.begin() ; object < objects.end(); object++ ){
    61 
    62                 /*Call echo on object: */
    63                 (*object)->Echo();
    64 
    65         }
    66         return;
    67 }
    68 
    69 void DataSet::DeepEcho(){
    70 
    71        
    72         vector<Object*>::iterator object;
    73 
    74         if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
    75 
    76         _printf_("DataSet echo: %i objects\n",objects.size());
    77 
    78         for ( object=objects.begin() ; object < objects.end(); object++ ){
    79 
    80                 /*Call deep echo on object: */
    81                 (*object)->DeepEcho();
    82 
    83         }
    84         return;
    85 }
    86 
    87 int DataSet::MarshallSize(){
    88 
    89         vector<Object*>::iterator object;
    90         int                      marshalled_dataset_size=0;
    91 
    92 
    93         for ( object=objects.begin() ; object < objects.end(); object++ ){
    94                 marshalled_dataset_size+= (*object)->MarshallSize();
    95         }
    96 
    97         marshalled_dataset_size+=sizeof(int); //objects size
    98         marshalled_dataset_size+=sizeof(int); //sorted size
    99         if(sorted){
    100                 marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids
    101                 marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets
    102         }
    103 
    104         return marshalled_dataset_size;
    105 }
    106 
     81/*}}}*/
     82
     83/*I/O*/
     84/*FUNCTION DataSet::Marshall{{{1*/
    10785char* DataSet::Marshall(){
    10886
    … …  
    150128        return marshalled_dataset;
    151129}
    152 
    153 int  DataSet::AddObject(Object* object){
    154 
    155         objects.push_back(object);
    156        
    157         return 1;
    158 }
    159 
    160 int  DataSet::DeleteObject(Object* object){
    161        
    162         vector<Object*>::iterator iterator;
    163 
    164         if(object){
    165                 iterator = find(objects.begin(), objects.end(),object);
    166                 delete *iterator;
    167                 objects.erase(iterator);
    168         }
    169 
    170 }
    171 
    172 int  DataSet::Size(void){
    173 
    174         return objects.size();
    175 }
    176 
     130/*}}}*/
     131/*FUNCTION DataSet::MarshallSize{{{1*/
     132int DataSet::MarshallSize(){
     133
     134        vector<Object*>::iterator object;
     135        int                      marshalled_dataset_size=0;
     136
     137
     138        for ( object=objects.begin() ; object < objects.end(); object++ ){
     139                marshalled_dataset_size+= (*object)->MarshallSize();
     140        }
     141
     142        marshalled_dataset_size+=sizeof(int); //objects size
     143        marshalled_dataset_size+=sizeof(int); //sorted size
     144        if(sorted){
     145                marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids
     146                marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets
     147        }
     148
     149        return marshalled_dataset_size;
     150}
     151/*}}}*/
     152/*FUNCTION DataSetDemarshall{{{1*/
    177153DataSet* DataSetDemarshall(char* marshalled_dataset){
    178154
    179155        int i;
    180        
     156
    181157        DataSet* dataset=NULL;
    182158        int      numobjects=0;
    … …  
    203179        }
    204180
    205         #ifdef _ISSM_DEBUG_
     181#ifdef _ISSM_DEBUG_
    206182        _printf_("Number of objects in dataset being demarshalled: %i\n",numobjects);
    207         #endif
     183#endif
    208184
    209185        for(i=0;i<numobjects;i++){
    210        
     186
    211187                /*get enum type of object: */
    212188                memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
    … …  
    310286
    311287}
    312 
     288/*}}}*/
     289
     290/*Specific methods*/
     291/*FUNCTION DataSet::AddObject{{{1*/
     292int  DataSet::AddObject(Object* object){
     293
     294        objects.push_back(object);
     295
     296        return 1;
     297}
     298/*}}}*/
     299/*FUNCTION DataSet::clear{{{1*/
     300void  DataSet::clear(){
     301
     302        vector<Object*>::iterator object;
     303
     304        for ( object=objects.begin() ; object < objects.end(); object++ ){
     305                delete (*object);
     306        }
     307        objects.clear();
     308}
     309/*}}}*/
     310/*FUNCTION DataSet::DeleteObject{{{1*/
     311int  DataSet::DeleteObject(Object* object){
     312
     313        vector<Object*>::iterator iterator;
     314
     315        if(object){
     316                iterator = find(objects.begin(), objects.end(),object);
     317                delete *iterator;
     318                objects.erase(iterator);
     319        }
     320
     321}
     322/*}}}*/
     323/*FUNCTION DataSet::DeepEcho{{{1*/
     324void DataSet::DeepEcho(){
     325
     326
     327        vector<Object*>::iterator object;
     328
     329        if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
     330
     331        _printf_("DataSet echo: %i objects\n",objects.size());
     332
     333        for ( object=objects.begin() ; object < objects.end(); object++ ){
     334
     335                /*Call deep echo on object: */
     336                (*object)->DeepEcho();
     337
     338        }
     339        return;
     340}
     341/*}}}*/
     342/*FUNCTION DataSet::Echo{{{1*/
     343void DataSet::Echo(){
     344
     345
     346        vector<Object*>::iterator object;
     347
     348        if(this==NULL)ISSMERROR(" trying to echo a NULL dataset");
     349
     350        _printf_("DataSet echo: %i objects\n",objects.size());
     351
     352        for ( object=objects.begin() ; object < objects.end(); object++ ){
     353
     354                /*Call echo on object: */
     355                (*object)->Echo();
     356
     357        }
     358        return;
     359}
     360/*}}}*/
     361/*FUNCTION DataSet::FindParam(double* pscalar, char* name){{{1*/
    313362int   DataSet::FindParam(double* pscalar, char* name){
    314363       
    … …  
    339388        return found;
    340389}
    341 
     390/*}}}*/
     391/*FUNCTION DataSet::FindParam(int* pinteger,char* name){{{1*/
    342392int   DataSet::FindParam(int* pinteger,char* name){
    343393       
    … …  
    369419        return found;
    370420}
    371 
     421/*}}}*/
     422/*FUNCTION DataSet::FindParam(char** pstring,char* name){{{1*/
    372423int   DataSet::FindParam(char** pstring,char* name){
    373424       
    … …  
    399450
    400451}
     452/*}}}*/
     453/*FUNCTION DataSet::FindParam(char*** pstringarray,int* pM,char* name){{{1*/
    401454int   DataSet::FindParam(char*** pstringarray,int* pM,char* name){
    402455       
    … …  
    429482
    430483}
     484/*}}}*/
     485/*FUNCTION DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){{{1*/
    431486int   DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){
    432487       
    … …  
    460515
    461516}
     517/*}}}*/
     518/*FUNCTION DataSet::FindParam(Vec* pvec,char* name){{{1*/
    462519int   DataSet::FindParam(Vec* pvec,char* name){
    463520       
    … …  
    489546
    490547}
     548/*}}}*/
     549/*FUNCTION DataSet::FindParamMat* pmat,char* name){{{1*/
    491550int   DataSet::FindParam(Mat* pmat,char* name){
    492551       
    … …  
    518577
    519578}
    520 
     579/*}}}*/
     580/*FUNCTION DataSet::FindParamObject{{{1*/
     581Object*   DataSet::FindParamObject(char* name){
     582
     583        /*Go through a dataset, and find a Param* object
     584         *which parameter name is "name" : */
     585
     586        vector<Object*>::iterator object;
     587        Param* param=NULL;
     588
     589        for ( object=objects.begin() ; object < objects.end(); object++ ){
     590
     591                /*Find param type objects: */
     592                if((*object)->Enum()==ParamEnum()){
     593
     594                        /*Ok, this object is a parameter, recover it and ask which name it has: */
     595                        param=(Param*)(*object);
     596
     597                        if (strcmp(param->GetParameterName(),name)==0){
     598                                /*Ok, this is the one! Return the object: */
     599                                return (*object);
     600                        }
     601                }
     602        }
     603        return NULL;
     604}
     605/*}}}*/
     606/*FUNCTION DataSet::FindResult(Vec* presult,char* name){{{1*/
    521607int   DataSet::FindResult(Vec* presult,char* name){
    522608
    … …  
    547633        return found;
    548634}
    549 
     635/*}}}*/
     636/*FUNCTION DataSet::FindResult(void* pvalue, char* name){{{1*/
    550637int   DataSet::FindResult(void* pvalue, char* name){
    551638
    … …  
    578665        return found;
    579666}
    580 
    581 
    582 Object*   DataSet::FindParamObject(char* name){
    583 
    584         /*Go through a dataset, and find a Param* object
    585          *which parameter name is "name" : */
    586        
    587         vector<Object*>::iterator object;
    588         Param* param=NULL;
    589 
    590         for ( object=objects.begin() ; object < objects.end(); object++ ){
    591 
    592                 /*Find param type objects: */
    593                 if((*object)->Enum()==ParamEnum()){
    594 
    595                         /*Ok, this object is a parameter, recover it and ask which name it has: */
    596                         param=(Param*)(*object);
    597 
    598                         if (strcmp(param->GetParameterName(),name)==0){
    599                                 /*Ok, this is the one! Return the object: */
    600                                 return (*object);
    601                         }
    602                 }
    603         }
    604         return NULL;
    605 }
    606 
    607 
     667/*}}}*/
     668/*FUNCTION DataSet::GetEnum(){{{1*/
     669int  DataSet::GetEnum(){
     670        return enum_type;
     671}
     672/*}}}*/
     673/*FUNCTION DataSet::GetEnum(int offset){{{1*/
     674int   DataSet::GetEnum(int offset){
     675
     676        return objects[offset]->Enum();
     677
     678}
     679/*}}}*/
     680/*FUNCTION DataSet::GetObjectByOffset{{{1*/
     681Object* DataSet::GetObjectByOffset(int offset){
     682
     683        return objects[offset];
     684
     685}
     686/*}}}*/
     687/*FUNCTION DataSet::GetObjectById{{{1*/
     688Object* DataSet::GetObjectById(int* poffset,int eid){
     689
     690        int id_offset;
     691        int offset;
     692        int i;
     693
     694        if(!sorted)ISSMERROR(" trying to binary search on a non-sorted dataset!");
     695
     696        /*Carry out a binary search on the sorted_ids: */
     697        if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){
     698                ISSMERROR(exprintf("%s%i"," could not find object with id ",eid));
     699        }
     700
     701        /*Convert  the id offset into sorted offset: */
     702        offset=id_offsets[id_offset];
     703
     704        /*Assign output pointers if requested:*/
     705        if (poffset)*poffset=offset;
     706
     707        /*Return object at offset position in objects :*/
     708        return objects[offset];
     709}
     710/*}}}*/
     711/*FUNCTION DataSet::NodeRank{{{1*/
    608712void   DataSet::NodeRank(int* ranks){
    609713
    610714        /*Go through a dataset, and find a Node* object. For this object,
    611715         * ask it to report its cpu rank: */
    612        
     716
    613717        int node_rank;
    614718        int id;
    615719        vector<Object*>::iterator object;
    616        
     720
    617721        for ( object=objects.begin() ; object < objects.end(); object++ ){
    618722
    619723                /*Call echo on object: */
    620724                if((*object)->Enum()==NodeEnum()){
    621                        
     725
    622726                        /*Ok, this object is a node, ask which rank it has: */
    623727                        node_rank=(*object)->MyRank();
    624                        
     728
    625729                        /*Which id does it have: */
    626730                        id=(*object)->GetId();
    … …  
    632736        return;
    633737}
    634 
    635 
     738/*}}}*/
     739/*FUNCTION DataSet::Presort{{{1*/
     740void DataSet::Presort(){
     741
     742        /*vector of objects is already sorted, just allocate the sorted ids and their
     743         * offsets:*/
     744        int i;
     745
     746        if(objects.size()){
     747                sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
     748                id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
     749                for(i=0;i<objects.size();i++){
     750                        id_offsets[i]=i;
     751                        sorted_ids[i]=objects[i]->GetId();
     752                }
     753        }
     754
     755        /*set sorted flag: */
     756        sorted=1;
     757}
     758/*}}}*/
     759/*FUNCTION DataSet::SetSorting{{{1*/
     760void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){
     761
     762        sorted=1;
     763        sorted_ids=in_sorted_ids;
     764        id_offsets=in_id_offsets;
     765}
     766/*}}}*/
     767/*FUNCTION DataSet::Size{{{1*/
     768int  DataSet::Size(void){
     769
     770        return objects.size();
     771}
     772/*}}}*/
     773/*FUNCTION DataSet::Sort{{{1*/
     774void DataSet::Sort(){
     775
     776        /*Only sort if we are not already sorted: */
     777        if(!sorted){
     778                ISSMERROR(" not implemented yet!");
     779        }
     780}
     781/*}}}*/
     782
     783/*Objects methods*/
     784/*FUNCTION DataSet::ComputePressure{{{1*/
     785void DataSet::ComputePressure(Vec p_g){
     786
     787        vector<Object*>::iterator object;
     788        Element* element=NULL;
     789
     790        for ( object=objects.begin() ; object < objects.end(); object++ ){
     791
     792                if(EnumIsElement((*object)->Enum())){
     793
     794                        element=(Element*)(*object);
     795                        element->ComputePressure(p_g);
     796                }
     797        }
     798
     799}
     800/*}}}*/
     801/*FUNCTION DataSet::Configure{{{1*/
     802void DataSet::Configure(DataSet* elements,DataSet* loads, DataSet* nodes, DataSet* materials,DataSet* parameters){
     803
     804        vector<Object*>::iterator object;
     805        Element* element=NULL;
     806        Load* load=NULL;
     807        Node* node=NULL;
     808        Numpar* numpar=NULL;
     809
     810        for ( object=objects.begin() ; object < objects.end(); object++ ){
     811
     812                if(EnumIsElement((*object)->Enum())){
     813
     814                        element=(Element*)(*object);
     815                        element->Configure(loads,nodes,materials,parameters);
     816                }
     817                if(EnumIsLoad((*object)->Enum())){
     818
     819                        load=(Load*)(*object);
     820
     821                        load->Configure(elements,nodes,materials);
     822                }
     823
     824                if((*object)->Enum()==NodeEnum()){
     825                        node=(Node*)(*object);
     826                        node->Configure(nodes);
     827                }
     828
     829                if((*object)->Enum()==NumparEnum()){
     830                        numpar=(Numpar*)(*object);
     831                        numpar->Configure(parameters);
     832                }
     833
     834        }
     835
     836}
     837/*}}}*/
     838/*FUNCTION DataSet::CostFunction{{{1*/
     839void  DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
     840
     841        double J=0;;
     842
     843        vector<Object*>::iterator object;
     844        Element* element=NULL;
     845
     846        for ( object=objects.begin() ; object < objects.end(); object++ ){
     847
     848                if(EnumIsElement((*object)->Enum())){
     849
     850                        element=(Element*)(*object);
     851                        J+=element->CostFunction(inputs,analysis_type,sub_analysis_type);
     852
     853                }
     854        }
     855
     856        /*Assign output pointers:*/
     857        *pJ=J;
     858
     859}
     860/*}}}*/
     861/*FUNCTION DataSet::CreateKMatrix{{{1*/
     862void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
     863
     864        vector<Object*>::iterator object;
     865        Element* element=NULL;
     866        Load* load=NULL;
     867
     868        for ( object=objects.begin() ; object < objects.end(); object++ ){
     869               
     870                if(EnumIsElement((*object)->Enum())){
     871
     872                        element=(Element*)(*object);
     873                        element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
     874                }
     875                if(EnumIsLoad((*object)->Enum())){
     876
     877                        load=(Load*)(*object);
     878                        load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
     879                }
     880        }
     881
     882}
     883/*}}}*/
     884/*FUNCTION DataSet::CreatePartitioningVector{{{1*/
     885void  DataSet::CreatePartitioningVector(Vec* ppartition,int numberofnodes,int numdofspernode){
     886
     887        /*output: */
     888        Vec partition=NULL;
     889        vector<Object*>::iterator object;
     890        Node* node=NULL;
     891
     892        /*Create partition vector: */
     893        partition=NewVec(numberofnodes);
     894
     895        /*Go through all nodes, and ask each node to plug its 1D doflist in
     896         * partition. The location where each node plugs its doflist into
     897         * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs).
     898         */
     899
     900        for ( object=objects.begin() ; object < objects.end(); object++ ){
     901
     902                /*Check this is a node: */
     903                if((*object)->Enum()==NodeEnum()){
     904
     905                        node=(Node*)(*object);
     906
     907                        /*Ok, this object is a node, ask it to plug values into partition: */
     908                        node->CreatePartition(partition);
     909
     910                }
     911        }
     912
     913        /*Assemble the petsc vector: */
     914        VecAssemblyBegin(partition);
     915        VecAssemblyEnd(partition);
     916
     917        /*Assign output pointers: */
     918        *ppartition=partition;
     919
     920        return;
     921}
     922/*}}}*/
     923/*FUNCTION DataSet::CreatePVector{{{1*/
     924void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
     925
     926        vector<Object*>::iterator object;
     927        Element* element=NULL;
     928        Load* load=NULL;
     929
     930        for ( object=objects.begin() ; object < objects.end(); object++ ){
     931
     932                if(EnumIsElement((*object)->Enum())){
     933
     934                        element=(Element*)(*object);
     935                        element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
     936                }
     937                if(EnumIsLoad((*object)->Enum())){
     938
     939                        load=(Load*)(*object);
     940                        load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
     941                }               
     942        }
     943
     944}
     945/*}}}*/
     946/*FUNCTION DataSet::DistributeDofs{{{1*/
    636947void  DataSet::DistributeDofs(int numberofnodes,int numdofspernode){
    637948
    … …  
    650961        extern int num_procs;
    651962        extern int my_rank;
    652        
     963
    653964        dofcount=0;
    654965        for ( object=objects.begin() ; object < objects.end(); object++ ){
    … …  
    658969
    659970                        node=(Node*)(*object);
    660                        
     971
    661972                        /*Ok, this object is a node, ask it to distribute dofs, and update the dofcount: */
    662973                        node->DistributeDofs(&dofcount,&dofcount1);
    663                        
     974
    664975                }
    665976        }
    … …  
    668979         *0. This means the dofs between all the cpus are not synchronized! We need to synchronize all
    669980         *dof on all cpus, and use those to update the dofs of every node: */
    670        
     981
    671982        alldofcount=(int*)xmalloc(num_procs*sizeof(int));
    672983        MPI_Gather(&dofcount,1,MPI_INT,alldofcount,1,MPI_INT,0,MPI_COMM_WORLD);
    … …  
    7001011
    7011012                        node=(Node*)(*object);
    702                        
     1013
    7031014                        /*Ok, this object is a node, ask it to update his dofs: */
    7041015                        node->UpdateDofs(dofcount,dofcount1);
    705                        
     1016
    7061017                }
    7071018        }
    … …  
    7201031
    7211032                        node=(Node*)(*object);
    722                        
     1033
    7231034                        /*Ok, let this object show its border dofs, if is is a border dof: */
    7241035                        node->ShowBorderDofs(borderdofs,borderdofs1);
    725                        
     1036
    7261037                }
    7271038        }
    … …  
    7361047
    7371048                        node=(Node*)(*object);
    738                        
     1049
    7391050                        /*Ok, this object is a node, ask it to update his dofs: */
    7401051                        node->UpdateBorderDofs(allborderdofs,allborderdofs1);
    741                        
     1052
    7421053                }
    7431054        }
    … …  
    7551066
    7561067}
    757 
    758                
    759 void  DataSet::CreatePartitioningVector(Vec* ppartition,int numberofnodes,int numdofspernode){
    760 
    761         /*output: */
    762         Vec partition=NULL;
     1068/*}}}*/
     1069/*FUNCTION DataSet::Du{{{1*/
     1070void  DataSet::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){
     1071
     1072
     1073        vector<Object*>::iterator object;
     1074        Element* element=NULL;
     1075
     1076        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1077
     1078                if(EnumIsElement((*object)->Enum())){
     1079
     1080                        element=(Element*)(*object);
     1081                        element->Du(du_g,inputs,analysis_type,sub_analysis_type);
     1082                }
     1083        }
     1084
     1085
     1086}
     1087/*}}}*/
     1088/*FUNCTION DataSet::FieldDepthAverageAtBase{{{1*/
     1089void  DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
     1090
    7631091        vector<Object*>::iterator object;
    7641092        Node* node=NULL;
    7651093
    766         /*Create partition vector: */
    767         partition=NewVec(numberofnodes);
    768 
    769         /*Go through all nodes, and ask each node to plug its 1D doflist in
    770          * partition. The location where each node plugs its doflist into
    771          * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs).
    772          */
    773        
    774         for ( object=objects.begin() ; object < objects.end(); object++ ){
    775 
    776                 /*Check this is a node: */
     1094        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1095
    7771096                if((*object)->Enum()==NodeEnum()){
    778                        
     1097
    7791098                        node=(Node*)(*object);
    780                        
    781                         /*Ok, this object is a node, ask it to plug values into partition: */
    782                         node->CreatePartition(partition);
    783                        
    784                 }
    785         }
    786 
    787         /*Assemble the petsc vector: */
    788         VecAssemblyBegin(partition);
    789         VecAssemblyEnd(partition);
    790 
    791         /*Assign output pointers: */
    792         *ppartition=partition;
    793 
    794         return;
    795 }
    796 
     1099                        node->FieldDepthAverageAtBase(field,field_serial,fieldname);
     1100                }
     1101        }
     1102
     1103}
     1104/*}}}*/
     1105/*FUNCTION DataSet::FieldExtrude{{{1*/
     1106void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
     1107
     1108        vector<Object*>::iterator object;
     1109        Penta* penta=NULL;
     1110        Node*  node=NULL;
     1111
     1112        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1113
     1114                if((*object)->Enum()==PentaEnum()){
     1115
     1116                        penta=(Penta*)(*object);
     1117                        penta->FieldExtrude(field,field_serial,field_name,collapse);
     1118
     1119                }
     1120                if((*object)->Enum()==NodeEnum()){
     1121
     1122                        node=(Node*)(*object);
     1123                        node->FieldExtrude(field,field_serial,field_name);
     1124
     1125                }
     1126        }
     1127
     1128}
     1129/*}}}*/
     1130/*FUNCTION DataSet::FlagClones{{{1*/
    7971131void  DataSet::FlagClones(int numberofnodes){
    7981132
    7991133        int i;
    8001134        extern int num_procs;
    801        
     1135
    8021136        int* ranks=NULL;
    8031137        int* minranks=NULL;
    … …  
    8051139        vector<Object*>::iterator object;
    8061140        Node* node=NULL;
    807        
     1141
    8081142        /*Allocate ranks: */
    8091143        ranks=(int*)xmalloc(numberofnodes*sizeof(int));
    8101144        minranks=(int*)xmalloc(numberofnodes*sizeof(int));
    811        
     1145
    8121146        for(i=0;i<numberofnodes;i++)ranks[i]=num_procs; //no cpu can have rank num_procs. This is the maximum limit.
    8131147
    … …  
    8151149        NodeRank(ranks);
    8161150
    817         #ifdef _ISSM_DEBUG_
     1151#ifdef _ISSM_DEBUG_
    8181152        for(i=0;i<numberofnodes;i++){
    8191153                _printf_("%i\n",ranks[i]);
    8201154        }
    821         #endif
     1155#endif
    8221156
    8231157        /*We need to take the minimum rank for each node, and every cpu needs to get that result. That way,
    … …  
    8271161        MPI_Allreduce ( (void*)ranks,(void*)minranks,numberofnodes,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
    8281162
    829         #ifdef _ISSM_DEBUG_
     1163#ifdef _ISSM_DEBUG_
    8301164        for(i=0;i<numberofnodes;i++){
    8311165                _printf_("%i\n",minranks[i]);
    8321166        }
    833         #endif
     1167#endif
    8341168
    8351169        /*Now go through all nodes, and use minranks to flag which nodes are cloned: */
    … …  
    8381172                /*Check this is a node: */
    8391173                if((*object)->Enum()==NodeEnum()){
    840                        
     1174
    8411175                        node=(Node*)(*object);
    842                        
     1176
    8431177                        /*Ok, this object is a node, ask it to plug values into partition: */
    8441178                        node->SetClone(minranks);
    845                        
     1179
    8461180                }
    8471181        }
    … …  
    8521186
    8531187}
     1188/*}}}*/
     1189/*FUNCTION DataSet::FlagNodeSets{{{1*/
     1190void DataSet::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s){
     1191
     1192        vector<Object*>::iterator object;
     1193        Node* node=NULL;
     1194
     1195        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1196
     1197                /*Check this is a single point constraint (spc): */
     1198                if((*object)->Enum()==NodeEnum()){
     1199
     1200                        node=(Node*)(*object);
     1201
     1202                        /*Plug set values intp our 4 set vectors: */
     1203                        node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s);
     1204
     1205                }
     1206        }
     1207
     1208        /*Assemble: */
     1209        VecAssemblyBegin(pv_g);
     1210        VecAssemblyEnd(pv_g);
     1211
     1212        VecAssemblyBegin(pv_m);
     1213        VecAssemblyEnd(pv_m);
     1214
     1215        VecAssemblyBegin(pv_n);
     1216        VecAssemblyEnd(pv_n);
     1217
     1218        VecAssemblyBegin(pv_f);
     1219        VecAssemblyEnd(pv_f);
     1220
     1221        VecAssemblyBegin(pv_s);
     1222        VecAssemblyEnd(pv_s);
     1223
     1224}
     1225/*}}}*/
     1226/*FUNCTION DataSet::Gradj{{{1*/
     1227void  DataSet::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
     1228
     1229
     1230        vector<Object*>::iterator object;
     1231        Element* element=NULL;
     1232
     1233        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1234
     1235                if(EnumIsElement((*object)->Enum())){
     1236
     1237                        element=(Element*)(*object);
     1238                        element->Gradj(grad_g,inputs,analysis_type,sub_analysis_type,control_type);
     1239                }
     1240        }
     1241
     1242
     1243}               
     1244/*}}}*/
     1245/*FUNCTION DataSet::MeltingIsPresent{{{1*/
     1246int   DataSet::MeltingIsPresent(){
     1247
     1248        int found=0;
     1249        int mpi_found=0;
     1250
     1251        vector<Object*>::iterator object;
     1252
     1253        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1254
     1255                if((*object)->Enum()==PengridEnum()){
     1256                        found=1;
     1257                        break;
     1258                }
     1259        }       
     1260       
     1261        #ifdef _PARALLEL_
     1262        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     1263        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     1264        found=mpi_found;
     1265        #endif
     1266
     1267        return found;
     1268}
     1269/*}}}*/
     1270/*FUNCTION DataSet::MeltingConstraints{{{1*/
     1271void  DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
     1272
     1273        /* generic object pointer: */
     1274        vector<Object*>::iterator object;
     1275        Pengrid* pengrid=NULL;
     1276
     1277        int unstable=0;
     1278        int num_unstable_constraints=0;
     1279        int converged=0;
     1280        int sum_num_unstable_constraints=0;
     1281
     1282        num_unstable_constraints=0;     
     1283
     1284        /*Enforce constraints: */
     1285        for ( object=objects.begin() ; object < objects.end(); object++ ){
    8541286               
    855 
     1287                if((*object)->Enum()==PengridEnum()){
     1288
     1289                        pengrid=(Pengrid*)(*object);
     1290
     1291                        pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type);
     1292
     1293                        num_unstable_constraints+=unstable;
     1294                }
     1295        }
     1296
     1297        #ifdef _PARALLEL_
     1298        MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     1299        MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);               
     1300        num_unstable_constraints=sum_num_unstable_constraints;
     1301        #endif
     1302
     1303        /*Have we converged? : */
     1304        if (num_unstable_constraints==0) converged=1;
     1305        else converged=0;
     1306
     1307        /*Assign output pointers: */
     1308        *pconverged=converged;
     1309        *pnum_unstable_constraints=num_unstable_constraints;
     1310}
     1311/*}}}*/
     1312/*FUNCTION DataSet::Misfit{{{1*/
     1313void  DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
     1314
     1315        double J=0;;
     1316
     1317        vector<Object*>::iterator object;
     1318        Element* element=NULL;
     1319
     1320        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1321
     1322                if(EnumIsElement((*object)->Enum())){
     1323
     1324                        element=(Element*)(*object);
     1325                        J+=element->Misfit(inputs,analysis_type,sub_analysis_type);
     1326
     1327                }
     1328        }
     1329
     1330        /*Assign output pointers:*/
     1331        *pJ=J;
     1332
     1333}
     1334/*}}}*/
     1335/*FUNCTION DataSet::NumberOfDofs{{{1*/
    8561336int   DataSet::NumberOfDofs(){
    8571337
    … …  
    8671347                /*Check this is a node: */
    8681348                if((*object)->Enum()==NodeEnum()){
    869                        
     1349
    8701350                        node=(Node*)(*object);
    871                        
     1351
    8721352                        /*Ok, this object is a node, ask it to plug values into partition: */
    8731353                        if (!node->IsClone()){
    8741354
    8751355                                numdofs+=node->GetNumberOfDofs();
    876                                
     1356
    8771357                        }
    8781358                }
    … …  
    8841364        return allnumdofs;
    8851365}
    886 
    887 void   DataSet::SetupSpcs(DataSet* nodes,Vec yg){
    888 
    889         vector<Object*>::iterator object;
    890         Spc* spc=NULL;
    891         Node* node=NULL;
    892 
    893         int nodeid;
    894         int dof;
    895         double value;
    896 
    897         for ( object=objects.begin() ; object < objects.end(); object++ ){
    898 
    899                 /*Check this is a single point constraint (spc): */
    900                 if((*object)->Enum()==SpcEnum()){
    901                        
    902                         spc=(Spc*)(*object);
    903                        
    904                         /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */
    905                         nodeid=spc->GetNodeId();
    906                         dof=spc->GetDof();
    907                         value=spc->GetValue();
    908 
    909                         /*Now, chase through nodes and find the corect node: */
    910                         node=(Node*)nodes->GetObjectById(NULL,nodeid);
    911 
    912                         /*Apply constraint: */
    913                         if(node){ //in case the spc is dealing with a node on another cpu
    914                                 node->ApplyConstraint(yg,dof,value);
    915                         }
    916 
    917                 }
    918         }
    919 
    920         /*Assemble yg: */
    921         VecAssemblyBegin(yg);
    922         VecAssemblyEnd(yg);
    923 }
    924                
    925 
    926                
     1366/*}}}*/
     1367/*FUNCTION DataSet::NumberOfRgbs{{{1*/
    9271368int   DataSet::NumberOfRgbs(){
    9281369
    … …  
    9431384        return count;
    9441385}
    945 
     1386/*}}}*/
     1387/*FUNCTION DataSet::OutputRifts{{{1*/
     1388void  DataSet::OutputRifts(Vec riftproperties){
     1389
     1390        vector<Object*>::iterator object;
     1391        Riftfront* riftfront=NULL;
     1392
     1393        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1394
     1395                if((*object)->Enum()==RiftfrontEnum()){
     1396
     1397                        riftfront=(Riftfront*)(*object);
     1398                        riftfront->OutputProperties(riftproperties);
     1399                }
     1400        }
     1401
     1402
     1403
     1404}
     1405/*}}}*/
     1406/*FUNCTION DataSet::PenaltyCreateKMatrix{{{1*/
     1407void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
     1408
     1409        vector<Object*>::iterator object;
     1410        Load* load=NULL;
     1411
     1412        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1413               
     1414                if(EnumIsLoad((*object)->Enum())){
     1415
     1416                        load=(Load*)(*object);
     1417                        load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
     1418                }
     1419        }
     1420
     1421}
     1422/*}}}*/
     1423/*FUNCTION DataSet::PenaltyCreatePVector{{{1*/
     1424void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
     1425
     1426        vector<Object*>::iterator object;
     1427        Load* load=NULL;
     1428
     1429        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1430
     1431                if(EnumIsLoad((*object)->Enum())){
     1432
     1433                        load=(Load*)(*object);
     1434                        load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
     1435                }               
     1436        }
     1437
     1438}
     1439/*}}}*/
     1440/*FUNCTION DataSet::RiftIsPresent{{{1*/
     1441int   DataSet::RiftIsPresent(){
     1442
     1443        int i;
     1444
     1445        vector<Object*>::iterator object;
     1446        Penpair* penpair=NULL;
     1447        int found=0;
     1448        int mpi_found=0;
     1449
     1450        /*go though loads, and figure out if one of the loads is a PenPair with numdof=2: */
     1451        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1452
     1453                if((*object)->Enum()==PenpairEnum()){
     1454
     1455                        penpair=(Penpair*)(*object);
     1456                }
     1457        }
     1458
     1459#ifdef _PARALLEL_
     1460        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     1461        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     1462        found=mpi_found;
     1463#endif
     1464
     1465        return found;
     1466}
     1467/*}}}*/
     1468/*FUNCTION DataSet::SetupMpcs{{{1*/
    9461469void DataSet::SetupMpcs(Mat Rmg,DataSet* nodes){
    9471470
    … …  
    9971520                                        node1->DofInMSet(dof-1);
    9981521                                }
    999                                
     1522
    10001523                                /*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the
    10011524                                 *same: */
    … …  
    10101533        }
    10111534}
    1012        
    1013 void DataSet::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s){
    1014 
    1015         vector<Object*>::iterator object;
     1535/*}}}*/
     1536/*FUNCTION DataSet::SetupSpcs{{{1*/
     1537void   DataSet::SetupSpcs(DataSet* nodes,Vec yg){
     1538
     1539        vector<Object*>::iterator object;
     1540        Spc* spc=NULL;
    10161541        Node* node=NULL;
    10171542
     1543        int nodeid;
     1544        int dof;
     1545        double value;
     1546
    10181547        for ( object=objects.begin() ; object < objects.end(); object++ ){
    10191548
    10201549                /*Check this is a single point constraint (spc): */
    1021                 if((*object)->Enum()==NodeEnum()){
    1022                        
    1023                         node=(Node*)(*object);
    1024 
    1025                         /*Plug set values intp our 4 set vectors: */
    1026                         node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s);
    1027 
    1028                 }
    1029         }
    1030 
    1031         /*Assemble: */
    1032         VecAssemblyBegin(pv_g);
    1033         VecAssemblyEnd(pv_g);
    1034 
    1035         VecAssemblyBegin(pv_m);
    1036         VecAssemblyEnd(pv_m);
    1037 
    1038         VecAssemblyBegin(pv_n);
    1039         VecAssemblyEnd(pv_n);
    1040 
    1041         VecAssemblyBegin(pv_f);
    1042         VecAssemblyEnd(pv_f);
    1043 
    1044         VecAssemblyBegin(pv_s);
    1045         VecAssemblyEnd(pv_s);
    1046 
    1047 }
    1048 void  DataSet::clear(){
    1049        
    1050         vector<Object*>::iterator object;
    1051        
    1052         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1053                 delete (*object);
    1054         }
    1055         objects.clear();
    1056 }
    1057 
    1058 
    1059 void DataSet::Configure(DataSet* elements,DataSet* loads, DataSet* nodes, DataSet* materials,DataSet* parameters){
     1550                if((*object)->Enum()==SpcEnum()){
     1551
     1552                        spc=(Spc*)(*object);
     1553
     1554                        /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */
     1555                        nodeid=spc->GetNodeId();
     1556                        dof=spc->GetDof();
     1557                        value=spc->GetValue();
     1558
     1559                        /*Now, chase through nodes and find the corect node: */
     1560                        node=(Node*)nodes->GetObjectById(NULL,nodeid);
     1561
     1562                        /*Apply constraint: */
     1563                        if(node){ //in case the spc is dealing with a node on another cpu
     1564                                node->ApplyConstraint(yg,dof,value);
     1565                        }
     1566
     1567                }
     1568        }
     1569
     1570        /*Assemble yg: */
     1571        VecAssemblyBegin(yg);
     1572        VecAssemblyEnd(yg);
     1573}
     1574/*}}}*/
     1575/*FUNCTION DataSet::SurfaceArea{{{1*/
     1576void  DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){
     1577
     1578        double S=0;;
    10601579
    10611580        vector<Object*>::iterator object;
    10621581        Element* element=NULL;
    1063         Load* load=NULL;
    1064         Node* node=NULL;
    1065         Numpar* numpar=NULL;
    1066 
    1067         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1068                
     1582
     1583        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1584
    10691585                if(EnumIsElement((*object)->Enum())){
    10701586
    10711587                        element=(Element*)(*object);
    1072                         element->Configure(loads,nodes,materials,parameters);
    1073                 }
    1074                 if(EnumIsLoad((*object)->Enum())){
    1075 
    1076                         load=(Load*)(*object);
    1077 
    1078                         load->Configure(elements,nodes,materials);
    1079                 }
    1080 
    1081                 if((*object)->Enum()==NodeEnum()){
    1082                         node=(Node*)(*object);
    1083                         node->Configure(nodes);
    1084                 }
    1085                
    1086                 if((*object)->Enum()==NumparEnum()){
    1087                         numpar=(Numpar*)(*object);
    1088                         numpar->Configure(parameters);
    1089                 }
    1090 
    1091         }
    1092 
    1093 }
    1094 
    1095 void DataSet::Presort(){
    1096        
    1097         /*vector of objects is already sorted, just allocate the sorted ids and their
    1098          * offsets:*/
    1099         int i;
    1100 
    1101         if(objects.size()){
    1102                 sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
    1103                 id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
    1104                 for(i=0;i<objects.size();i++){
    1105                         id_offsets[i]=i;
    1106                         sorted_ids[i]=objects[i]->GetId();
    1107                 }
    1108         }
    1109 
    1110         /*set sorted flag: */
    1111         sorted=1;
    1112 }
    1113 
    1114 void DataSet::Sort(){
    1115 
    1116         /*Only sort if we are not already sorted: */
    1117         if(!sorted){
    1118                 ISSMERROR(" not implemented yet!");
    1119         }
    1120 }
    1121 
    1122 
    1123 Object* DataSet::GetObjectByOffset(int offset){
    1124 
    1125         return objects[offset];
    1126 
    1127 }
    1128 
    1129 Object* DataSet::GetObjectById(int* poffset,int eid){
    1130 
    1131         int id_offset;
    1132         int offset;
    1133         int i;
    1134 
    1135         if(!sorted)ISSMERROR(" trying to binary search on a non-sorted dataset!");
    1136 
    1137         /*Carry out a binary search on the sorted_ids: */
    1138         if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){
    1139                 ISSMERROR(exprintf("%s%i"," could not find object with id ",eid));
    1140         }
    1141 
    1142         /*Convert  the id offset into sorted offset: */
    1143         offset=id_offsets[id_offset];
    1144 
    1145         /*Assign output pointers if requested:*/
    1146         if (poffset)*poffset=offset;
    1147 
    1148         /*Return object at offset position in objects :*/
    1149         return objects[offset];
    1150 }
    1151 
    1152 void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){
    1153 
    1154         sorted=1;
    1155         sorted_ids=in_sorted_ids;
    1156         id_offsets=in_id_offsets;
    1157 }
    1158 
    1159 void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    1160 
    1161         vector<Object*>::iterator object;
    1162         Element* element=NULL;
    1163         Load* load=NULL;
    1164 
    1165         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1166                
    1167                 if(EnumIsElement((*object)->Enum())){
    1168 
    1169                         element=(Element*)(*object);
    1170                         element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
    1171                 }
    1172                 if(EnumIsLoad((*object)->Enum())){
    1173 
    1174                         load=(Load*)(*object);
    1175                         load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
    1176                 }
    1177         }
    1178 
    1179 }
    1180 
    1181 void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    1182 
    1183         vector<Object*>::iterator object;
    1184         Element* element=NULL;
    1185         Load* load=NULL;
    1186 
    1187         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1188 
    1189                 if(EnumIsElement((*object)->Enum())){
    1190 
    1191                         element=(Element*)(*object);
    1192                         element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    1193                 }
    1194                 if(EnumIsLoad((*object)->Enum())){
    1195 
    1196                         load=(Load*)(*object);
    1197                         load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
    1198                 }               
    1199         }
    1200 
    1201 }
    1202 
     1588                        S+=element->SurfaceArea(inputs,analysis_type,sub_analysis_type);
     1589
     1590                }
     1591        }
     1592
     1593        /*Assign output pointers:*/
     1594        *pS=S;
     1595
     1596}
     1597/*}}}*/
     1598/*FUNCTION DataSet::UpdateFromInputs{{{1*/
    12031599void  DataSet::UpdateFromInputs(void* inputs){
    12041600
    … …  
    12331629
    12341630}
    1235 
    1236 void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    1237 
    1238         vector<Object*>::iterator object;
    1239         Load* load=NULL;
    1240 
    1241         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1242                
    1243                 if(EnumIsLoad((*object)->Enum())){
    1244 
    1245                         load=(Load*)(*object);
    1246                         load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
    1247                 }
    1248         }
    1249 
    1250 }
    1251 
    1252 void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    1253 
    1254         vector<Object*>::iterator object;
    1255         Load* load=NULL;
    1256 
    1257         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1258 
    1259                 if(EnumIsLoad((*object)->Enum())){
    1260 
    1261                         load=(Load*)(*object);
    1262                         load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
    1263                 }               
    1264         }
    1265 
    1266 }
    1267 
    1268 int   DataSet::RiftIsPresent(){
    1269 
    1270         int i;
    1271        
    1272         vector<Object*>::iterator object;
    1273         Penpair* penpair=NULL;
    1274         int found=0;
    1275         int mpi_found=0;
    1276 
    1277         /*go though loads, and figure out if one of the loads is a PenPair with numdof=2: */
    1278         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1279 
    1280                 if((*object)->Enum()==PenpairEnum()){
    1281 
    1282                         penpair=(Penpair*)(*object);
    1283                 }
    1284         }
    1285 
    1286         #ifdef _PARALLEL_
    1287         MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
    1288         MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
    1289         found=mpi_found;
    1290         #endif
    1291 
    1292         return found;
    1293 }
    1294 
    1295 int   DataSet::MeltingIsPresent(){
    1296 
    1297         int found=0;
    1298         int mpi_found=0;
    1299 
    1300         vector<Object*>::iterator object;
    1301 
    1302         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1303 
    1304                 if((*object)->Enum()==PengridEnum()){
    1305                         found=1;
    1306                         break;
    1307                 }
    1308         }       
    1309        
    1310         #ifdef _PARALLEL_
    1311         MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
    1312         MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
    1313         found=mpi_found;
    1314         #endif
    1315 
    1316         return found;
    1317 }
    1318 
    1319 void  DataSet::MeltingConstraints(int* pconverged, int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
    1320 
    1321         /* generic object pointer: */
    1322         vector<Object*>::iterator object;
    1323         Pengrid* pengrid=NULL;
    1324 
    1325         int unstable=0;
    1326         int num_unstable_constraints=0;
    1327         int converged=0;
    1328         int sum_num_unstable_constraints=0;
    1329 
    1330         num_unstable_constraints=0;     
    1331 
    1332         /*Enforce constraints: */
    1333         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1334                
    1335                 if((*object)->Enum()==PengridEnum()){
    1336 
    1337                         pengrid=(Pengrid*)(*object);
    1338 
    1339                         pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type);
    1340 
    1341                         num_unstable_constraints+=unstable;
    1342                 }
    1343         }
    1344 
    1345         #ifdef _PARALLEL_
    1346         MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
    1347         MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);               
    1348         num_unstable_constraints=sum_num_unstable_constraints;
    1349         #endif
    1350 
    1351         /*Have we converged? : */
    1352         if (num_unstable_constraints==0) converged=1;
    1353         else converged=0;
    1354 
    1355         /*Assign output pointers: */
    1356         *pconverged=converged;
    1357         *pnum_unstable_constraints=num_unstable_constraints;
    1358 }
    1359 
    1360 DataSet*   DataSet::Copy(void){
    1361 
    1362 
    1363         DataSet* copy=NULL;
    1364         vector<Object*>::iterator object;
    1365         Object* object_copy=NULL;
    1366 
    1367         copy=new DataSet(enum_type);
    1368 
    1369         copy->sorted=sorted;
    1370         copy->presorted=presorted;
    1371         if(sorted_ids){
    1372                 copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
    1373                 memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int));
    1374         }
    1375         if(id_offsets){
    1376                 copy->id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
    1377                 memcpy(copy->id_offsets,id_offsets,objects.size()*sizeof(int));
    1378         }
    1379 
    1380         /*Now we need to deep copy the objects: */
    1381         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1382 
    1383                 /*Call echo on object: */
    1384                 object_copy = (*object)->copy();
    1385                 copy->AddObject(object_copy);
    1386         }
    1387         return copy;
    1388 }
    1389 
    1390 
    1391 void  DataSet::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){
    1392 
    1393 
    1394         vector<Object*>::iterator object;
    1395         Element* element=NULL;
    1396 
    1397         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1398                
    1399                 if(EnumIsElement((*object)->Enum())){
    1400 
    1401                         element=(Element*)(*object);
    1402                         element->Du(du_g,inputs,analysis_type,sub_analysis_type);
    1403                 }
    1404         }
    1405 
    1406 
    1407 }
    1408 
    1409 void  DataSet::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    1410 
    1411 
    1412         vector<Object*>::iterator object;
    1413         Element* element=NULL;
    1414 
    1415         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1416                
    1417                 if(EnumIsElement((*object)->Enum())){
    1418 
    1419                         element=(Element*)(*object);
    1420                         element->Gradj(grad_g,inputs,analysis_type,sub_analysis_type,control_type);
    1421                 }
    1422         }
    1423 
    1424 
    1425 }               
    1426                
    1427 void  DataSet::Misfit(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
    1428 
    1429         double J=0;;
    1430        
    1431         vector<Object*>::iterator object;
    1432         Element* element=NULL;
    1433 
    1434         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1435                
    1436                 if(EnumIsElement((*object)->Enum())){
    1437 
    1438                         element=(Element*)(*object);
    1439                         J+=element->Misfit(inputs,analysis_type,sub_analysis_type);
    1440 
    1441                 }
    1442         }
    1443 
    1444         /*Assign output pointers:*/
    1445         *pJ=J;
    1446 
    1447 }
    1448 
    1449 void  DataSet::CostFunction(double* pJ,void* inputs,int analysis_type,int sub_analysis_type){
    1450 
    1451         double J=0;;
    1452 
    1453         vector<Object*>::iterator object;
    1454         Element* element=NULL;
    1455 
    1456         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1457 
    1458                 if(EnumIsElement((*object)->Enum())){
    1459 
    1460                         element=(Element*)(*object);
    1461                         J+=element->CostFunction(inputs,analysis_type,sub_analysis_type);
    1462 
    1463                 }
    1464         }
    1465 
    1466         /*Assign output pointers:*/
    1467         *pJ=J;
    1468 
    1469 }
    1470 
    1471 void  DataSet::SurfaceArea(double* pS,void* inputs,int analysis_type,int sub_analysis_type){
    1472 
    1473         double S=0;;
    1474 
    1475         vector<Object*>::iterator object;
    1476         Element* element=NULL;
    1477 
    1478         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1479 
    1480                 if(EnumIsElement((*object)->Enum())){
    1481 
    1482                         element=(Element*)(*object);
    1483                         S+=element->SurfaceArea(inputs,analysis_type,sub_analysis_type);
    1484 
    1485                 }
    1486         }
    1487 
    1488         /*Assign output pointers:*/
    1489         *pS=S;
    1490 
    1491 }
    1492 
    1493 void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
    1494 
    1495         vector<Object*>::iterator object;
    1496         Penta* penta=NULL;
    1497         Node*  node=NULL;
    1498 
    1499         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1500                
    1501                 if((*object)->Enum()==PentaEnum()){
    1502 
    1503                         penta=(Penta*)(*object);
    1504                         penta->FieldExtrude(field,field_serial,field_name,collapse);
    1505 
    1506                 }
    1507                 if((*object)->Enum()==NodeEnum()){
    1508 
    1509                         node=(Node*)(*object);
    1510                         node->FieldExtrude(field,field_serial,field_name);
    1511 
    1512                 }
    1513         }
    1514 
    1515 }
    1516 
    1517 void  DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
    1518 
    1519         vector<Object*>::iterator object;
    1520         Node* node=NULL;
    1521 
    1522         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1523 
    1524                 if((*object)->Enum()==NodeEnum()){
    1525 
    1526                         node=(Node*)(*object);
    1527                         node->FieldDepthAverageAtBase(field,field_serial,fieldname);
    1528                 }
    1529         }
    1530 
    1531 }
    1532 
    1533 void DataSet::ComputePressure(Vec p_g){
    1534 
    1535         vector<Object*>::iterator object;
    1536         Element* element=NULL;
    1537 
    1538         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1539                
    1540                 if(EnumIsElement((*object)->Enum())){
    1541 
    1542                         element=(Element*)(*object);
    1543                         element->ComputePressure(p_g);
    1544                 }
    1545         }
    1546 
    1547 }
    1548 
    1549 
     1631/*}}}*/
     1632/*FUNCTION DataSet::UpdateNodePositions{{{1*/
    15501633void  DataSet::UpdateNodePositions(double* thickness,double* bed){
    15511634
    … …  
    15621645        }
    15631646}
    1564                
    1565 
    1566 int   DataSet::GetEnum(int offset){
    1567 
    1568         return objects[offset]->Enum();
    1569 
    1570 }
    1571 
    1572 void  DataSet::OutputRifts(Vec riftproperties){
    1573 
    1574         vector<Object*>::iterator object;
    1575         Riftfront* riftfront=NULL;
    1576 
    1577         for ( object=objects.begin() ; object < objects.end(); object++ ){
    1578 
    1579                 if((*object)->Enum()==RiftfrontEnum()){
    1580 
    1581                         riftfront=(Riftfront*)(*object);
    1582                         riftfront->OutputProperties(riftproperties);
    1583                 }
    1584         }
    1585 
    1586 
    1587 
    1588 }
     1647/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.