Changeset 1628 for issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

Timestamp:

08/10/09 17:16:00 (16 years ago)

Author:

Eric.Larour

Message:

Added rifts formulation, in diagnostic 2d

File:

• issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (modified) (6 diffs)

issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

@@ -16 +16 @@
 void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){
 
-        int i,j,counter;
+        int i;
         int element;
 
@@ -24 +24 @@
         DataSet*    loads    = NULL;
         Icefront*   icefront = NULL;
-        Penpair*    penpair  = NULL;
+        Riftfront*  riftfront= NULL;
 
         int segment_width;
@@ -30 +30 @@
         int i1,i2,i3,i4;
         int i0,range;
-        int grid1,grid2;
-
-        /*rift penpair: */
-        double  normal[2];
-        double  length;
-        double  friction;
-        int     fill;
+
                 
         /*icefront intermediary data: */
@@ -48 +42 @@
         double  icefront_b[MAX_ICEFRONT_GRIDS];
 
-        /*penpair intermediary data: */
-        int penpair_id;
-        int penpair_node_ids[2];
-        double penpair_penalty_offset;
-        int penpair_numdofs;
-        int penpair_dof;
-        int penpair_penalty_lock;
-        int penpair_element_ids[2];
-        double penpair_friction;
-        int penpair_fill;
-        double penpair_normal[2];
-        double penpair_length;
-
-        /*Rifts:*/
-        int* riftsnumpenaltypairs=NULL;
-        double** riftspenaltypairs=NULL;
-        int* riftsfill=NULL;
-        int* riftsfriction=NULL;
-        double* riftpenaltypairs=NULL;
+        /*rifts: */
+        char riftfront_type[RIFTFRONTSTRING];
+        int  riftfront_id;
+        int  riftfront_node_ids[2];
+        int  riftfront_mparid;
+        double riftfront_h[2];
+        double riftfront_b[2];
+        double riftfront_s[2];
+        double riftfront_normal[2];
+        double riftfront_length;
+        int    riftfront_fill;
+        double riftfront_friction;
+        bool   riftfront_shelf;
+        double riftfront_penalty_offset;
+        int riftfront_penalty_lock;
+        bool riftfront_active;
+        int  riftfront_counter;
+        bool riftfront_prestable;
         int el1,el2;
+        int grid1,grid2;
+        double normal[2];
+        double length;
+        int    fill;
+        double friction;
         
         /*Create loads: */
@@ -153 +151 @@
         xfree((void**)&model->bed);
 
-        counter=0;
-
         /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
         /*First fetch rifts: */
+        ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
+        ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
+        ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
+        ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
+        ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
+
         if(model->numrifts){
-
                 for(i=0;i<model->numrifts;i++){
-                        riftpenaltypairs=model->riftspenaltypairs[i];
-                        for(j=0;j<model->riftsnumpenaltypairs[i];j++){
-
-                                el1=(int)*(riftpenaltypairs+7*j+2);
-                                #ifdef _PARALLEL_
-                                if (model->epart[el1-1]!=my_rank){
-                                        /*This load does not belong to this cluster node, as it references an element 
-                                         *that does not belong to this node's partition. Drop this 'j':*/
-                                        continue;
-                                }
-                                #endif
-
-                                /*Ok, retrieve all the data needed to add a penalty between the two grids: */
-                                el2=(int)*(riftpenaltypairs+7*j+3); 
-                                grid1=(int)*(riftpenaltypairs+7*j+0); 
-                                grid2=(int)*(riftpenaltypairs+7*j+1);
-                                normal[0]=*(riftpenaltypairs+7*j+4);
-                                normal[1]=*(riftpenaltypairs+7*j+5);
-                                length=*(riftpenaltypairs+7*j+6);
-                                friction=model->riftsfriction[i];
-                                fill=model->riftsfill[i];
-
-                                penpair_id=counter+1; //matlab indexing
-                                penpair_node_ids[0]=grid1;
-                                penpair_node_ids[1]=grid2;
-                                penpair_element_ids[0]=el1;
-                                penpair_element_ids[1]=el2;
-                                penpair_penalty_offset=model->penalty_offset;
-                                penpair_penalty_lock=model->penalty_lock;
-                                penpair_numdofs=2;
-                                penpair_normal[0]=normal[0];
-                                penpair_normal[1]=normal[1];
-                                penpair_length=length;
-                                penpair_friction=friction;
-                                penpair_fill=fill;
-
-                                penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
-                                                penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
                                 
-                                loads->AddObject(penpair);
-
-                                counter++;
+                        el1=(int)*(model->riftinfo+9*i+2);
+                        #ifdef _PARALLEL_
+                        if (model->epart[el1-1]!=my_rank){
+                                /*This load does not belong to this cluster node, as it references an element 
+                                 *that does not belong to this node's partition. Drop this 'j':*/
+                                continue;
                         }
+                        #endif
+
+                        /*Ok, retrieve all the data needed to add a penalty between the two grids: */
+                        el2=(int)*(model->riftinfo+9*i+3); 
+
+                        grid1=(int)*(model->riftinfo+9*i+0); 
+                        grid2=(int)*(model->riftinfo+9*i+1);
+                        
+                        normal[0]=*(model->riftinfo+9*i+4);
+                        normal[1]=*(model->riftinfo+9*i+5);
+                        length=*(model->riftinfo+9*i+6);
+                        
+                        fill = (int)*(model->riftinfo+9*i+7);
+                        friction=*(model->riftinfo+9*i+8);
+        
+                        strcpy(riftfront_type,"2d");
+                        riftfront_id=i+1; //matlab indexing
+                        riftfront_node_ids[0]=grid1;
+                        riftfront_node_ids[1]=grid2;
+                        riftfront_mparid=model->numberofelements+1; //matlab indexing
+
+                        riftfront_h[0]=model->thickness[grid1-1];
+                        riftfront_h[1]=model->thickness[grid2-1];
+
+                        riftfront_b[0]=model->bed[grid1-1];
+                        riftfront_b[1]=model->bed[grid2-1];
+
+                        riftfront_s[0]=model->surface[grid1-1];
+                        riftfront_s[1]=model->surface[grid2-1];
+
+                        riftfront_normal[0]=normal[0];
+                        riftfront_normal[1]=normal[1];
+                        riftfront_length=length;
+                        
+                        riftfront_fill=fill;
+                        riftfront_friction=friction;
+                        riftfront_shelf=(bool)model->gridoniceshelf[grid1-1];
+
+                        riftfront_penalty_offset=model->penalty_offset;
+                        riftfront_penalty_lock=model->penalty_lock;
+
+                        riftfront_active=0;
+                        riftfront_counter=0;
+                        riftfront_prestable=0;
+                        
+                        riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
+
+                        loads->AddObject(riftfront);
                 }
         }
+
+        /*free ressources: */
+        xfree((void**)&model->riftinfo);
+        xfree((void**)&model->thickness);
+        xfree((void**)&model->bed);
+        xfree((void**)&model->surface);
+        xfree((void**)&model->gridoniceshelf);
+
+        cleanup_and_return:
 
         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
@@ -210 +233 @@
         loads->Presort();
 
-
-        /*Free ressources:*/
-        xfree((void**)&riftsnumpenaltypairs); 
-        for(i=0;i<model->numrifts;i++){
-                double* temp=riftspenaltypairs[i];
-                xfree((void**)&temp);
-        }
-        xfree((void**)&riftspenaltypairs);
-        xfree((void**)&riftsfill);
-        xfree((void**)&riftsfriction);
-        
-        cleanup_and_return:
-
         /*Assign output pointer: */
         *ploads=loads;