Changeset 1196

Timestamp:

07/01/09 11:11:40 (16 years ago)

Author:

Mathieu Morlighem

Message:

Added Control method macayeal 3d

Location:

issm/trunk/src/c

Files:

• objects/Penta.cpp (modified) (5 diffs)
• parallel/GradJCompute.cpp (modified) (3 diffs)

issm/trunk/src/c/objects/Penta.cpp

@@ -1222 +1222 @@
 void  Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){
         
-        
         Tria* tria=NULL;
         
@@ -1241 +1240 @@
 #define __FUNCT__ "Penta::GradjB"
 void  Penta::GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){
-        throw ErrorException(__FUNCT__," not supported yet!");
+
+        Tria* tria=NULL;
+
+        if (collapse){
+                /*Bail out element if collapsed (2d) and not on bed*/
+                if (!onbed) return;
+
+                /*This element should be collapsed into a tria element at its base. Create this tria element, 
+                 * and compute gardj*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
+                tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type);
+                delete tria;
+                return;
+        }
+        else{
+                throw ErrorException(__FUNCT__," not supported yet!");
+        }
 }
         
@@ -1888 +1903 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__ "Penta::CreateKMatrix"
+#define __FUNCT__ "Penta::GetParameterValue"
 void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
         
@@ -2011 +2026 @@
                         }
                 } //if (strcmp(field_name,"velocity")==0)
+                else if (strcmp(field_name,"gradj")==0){
+
+                        /* node data: */
+                        int          dof1;
+                        double       fieldel;
+
+                        /*this penta is a collapsed macayeal. For each node on the base of this penta, 
+                         * we grab the field. Once we know the field, we follow the upper nodes, 
+                         * inserting the same field value into field, until we reach the surface: */
+                        for(i=0;i<3;i++){
+
+                                node=nodes[i]; //base nodes
+                                dof1=node->GetDofList1();
+
+                                /*get field for this base node: */
+                                fieldel=field_serial[dof1];
+
+                                //go throfieldn all nodes which sit on top of this node, until we reach the surface, 
+                                //and plfield  field in field
+                                for(;;){
+
+                                        dof1=node->GetDofList1();
+                                        VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
+
+                                        if (node->IsOnSurface())break;
+                                        /*get next node: */
+                                        node=node->GetUpperNode();
+                                }
+                        }
+                }
                 else if ( 
                                 (strcmp(field_name,"thickness")==0) ||
@@ -2016 +2061 @@
                                 (strcmp(field_name,"bed")==0)  ||
                                 (strcmp(field_name,"slopex")==0)  ||
-                                (strcmp(field_name,"slopey")==0) 
+                                (strcmp(field_name,"slopey")==0)
                                 ){
 

issm/trunk/src/c/parallel/GradJCompute.cpp

@@ -41 +41 @@
         Vec grad_g=NULL;
 
+        /*flags: */
+        int debug=0;
+        int dim=-1;
         
         /*some parameters:*/
@@ -48 +51 @@
         femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
         femmodel->parameters->FindParam((void*)&control_type,"control_type");
+        femmodel->parameters->FindParam((void*)&debug,"debug");
+        femmodel->parameters->FindParam((void*)&dim,"dim");
 
         _printf_("%s\n","      recover solution for this stiffness and right hand side:");
@@ -75 +80 @@
         Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, 
                                 inputs,analysis_type,sub_analysis_type,control_type);
+
+        if (dim==3){
+
+                _printf_("%s\n","      extruding gradient...");
+                FieldExtrudex( grad_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,"gradj",1);
+        }
         
         /*Free ressources:*/