Changeset 848 for issm/trunk/src/c/objects/Node.cpp

Timestamp:

06/08/09 15:52:02 (15 years ago)

Author:

Eric.Larour

Message:

Added FieldExtrudex, FieldDepthAveragex, UpdateNodePositionsx modules

File:

• issm/trunk/src/c/objects/Node.cpp (modified) (3 diffs)

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

@@ -550 +550 @@
                 
 #undef __FUNCT__ 
-#define __FUNCT__ "Node::VelocityDepthAverageAtBase"
-void  Node::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
+#define __FUNCT__ "Node::FieldDepthAverageAtBase"
+void  Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){
 
         /* node data: */
@@ -560 +560 @@
         Node* node=NULL;
         Node* upper_node=NULL;
-        double velocity2[2];
-        double velocity1[2];
-        double velocity_average[2];
-        double sum[2];
         double z1,z2,dz;
         double thickness;
 
+        /*Are we on the base, not on the surface, and not on a clone node?:*/
         
         if(onbed==1 & clone==0 &onsurface==0){
@@ -572 +569 @@
                 doflist1=GetDofList1();
 
-                /*this node is on the bed. We are going to:
-                 * follow the upper nodes until we reach the surface. At each upper node, we'll grab the 
-                 * velocity for this node, and add it to ug. We'll finish by 
-                 * we grab the velocity. Once we know the velocity, we follow the upper nodes, 
-                 * inserting the same velocity value into ug, until we reach the surface: */
-                sum[0]=0;
-                sum[1]=0;
-                thickness=0;
-
-                /*get dofs for this base node velocity: we know there are two dofs in ug_serial */
-                dofx=2*doflist1;
-                dofy=2*doflist1+1;
-
-                node=this;
-                for(;;){
-
-                        if (node->IsOnSurface())break;
-
-                        doflist1=node->GetDofList1();
+                /*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, 
+                 * we'll grab the * field for this node, and add it to overall field: */
+
+                if(strcmp(fieldname,"velocity")==0){
+
+                        /*field is a velocity, 2 dofs per node: */
+                        double velocity2[2];
+                        double velocity1[2];
+                        double velocity_average[2];
+                        double sum[2];
+
+                        sum[0]=0;
+                        sum[1]=0;
+                        thickness=0;
+
+                        /*get dofs for this base node velocity: we know there are two dofs in field_serial */
+                        dofx=2*doflist1;
+                        dofy=2*doflist1+1;
+
+                        node=this;
+                        for(;;){
+
+                                if (node->IsOnSurface())break;
+
+                                doflist1=node->GetDofList1();
+                                
+                                velocity1[0]=field_serial[2*doflist1];
+                                velocity1[1]=field_serial[2*doflist1+1];
+                                z1=node->GetZ();
+
+                                upper_node=node->GetUpperNode();
+                                doflist1=upper_node->GetDofList1();
                         
-                        velocity1[0]=ug_serial[2*doflist1];
-                        velocity1[1]=ug_serial[2*doflist1+1];
-                        z1=node->GetZ();
-
-                        upper_node=node->GetUpperNode();
-                        doflist1=upper_node->GetDofList1();
-                
-                        velocity2[0]=ug_serial[2*doflist1];
-                        velocity2[1]=ug_serial[2*doflist1+1];
-                        z2=upper_node->GetZ();
-
-                        dz=(z2-z1);
-                        thickness+=dz;
-                        velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
-                        velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
-
-                        sum[0]+=velocity_average[0]*dz;
-                        sum[1]+=velocity_average[1]*dz;
-
-                        /* get next node: */
-                        node=node->GetUpperNode();
+                                velocity2[0]=field_serial[2*doflist1];
+                                velocity2[1]=field_serial[2*doflist1+1];
+                                z2=upper_node->GetZ();
+
+                                dz=(z2-z1);
+                                thickness+=dz;
+                                velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
+                                velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
+
+                                sum[0]+=velocity_average[0]*dz;
+                                sum[1]+=velocity_average[1]*dz;
+
+                                /* get next node: */
+                                node=node->GetUpperNode();
+                        }
+
+                        sum[0]=sum[0]/thickness;
+                        sum[1]=sum[1]/thickness;
+
+                        /* Plfield velocity_average*deltaH/H into base of field: */
+                        VecSetValues(field,1,&dofx,&sum[0],INSERT_VALUES);
+                        VecSetValues(field,1,&dofy,&sum[1],INSERT_VALUES);
                 }
-
-                sum[0]=sum[0]/thickness;
-                sum[1]=sum[1]/thickness;
-
-                /* Plug velocity_average*deltaH/H into base of ug: */
-                VecSetValues(ug,1,&dofx,&sum[0],INSERT_VALUES);
-                VecSetValues(ug,1,&dofy,&sum[1],INSERT_VALUES);
-        }
-}
+                else{
+                        /*field is regular, 1 dof per node: */
+                        double field2;
+                        double field1;
+                        double field_average;
+                        double sum;
+
+                        sum=0;
+                        thickness=0;
+
+                        /*get dofs for this base node velocity: we know there are two dofs in field_serial */
+                        dofx=doflist1;
+
+                        node=this;
+                        for(;;){
+
+                                if (node->IsOnSurface())break;
+
+                                doflist1=node->GetDofList1();
+                                
+                                field1=field_serial[doflist1];
+                                z1=node->GetZ();
+
+                                upper_node=node->GetUpperNode();
+                                doflist1=upper_node->GetDofList1();
+                        
+                                field2=field_serial[2*doflist1];
+                                z2=upper_node->GetZ();
+
+                                dz=(z2-z1);
+                                thickness+=dz;
+                                field_average=(field1+field2)/2.0;
+
+                                sum+=field_average*dz;
+
+                                /* get next node: */
+                                node=node->GetUpperNode();
+                        }
+
+                        sum=sum/thickness;
+
+                        /* Plug field_average*deltH/H into base of field: */
+                        VecSetValues(field,1,&dofx,&sum,INSERT_VALUES);
+                }
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Node::FieldExtrude"
+void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
+                
+        /* node data: */
+        int   numberofdofspernode;
+        Node* node=NULL;
+        int   i;
+
+        /*Is this node on bed? :*/
+        if (onbed){
+
+                if (strcmp(field_name,"velocity")==0){
+
+                        /* node data: */
+                        const int    numdof=2;
+                        int          doflist[numdof];
+                        int          nodedofs[2];
+                        double       fieldnode[2];
+
+
+                        GetDofList(&doflist[0],&numberofdofspernode);
+
+                        //initilaize node
+                        node=this;
+                
+                        /*get field for this base node: */
+                        fieldnode[0]=field_serial[doflist[0]];
+                        fieldnode[1]=field_serial[doflist[1]];
+
+                        //go through all nodes which sit on top of this node, until we reach the surface, 
+                        //and plfield  field in field
+                        for(;;){
+
+                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
+                                VecSetValues(field,1,&nodedofs[0],&fieldnode[0],INSERT_VALUES);
+                                VecSetValues(field,1,&nodedofs[1],&fieldnode[1],INSERT_VALUES);
+
+                                if (node->IsOnSurface())break;
+                                /*get next node: */
+                                node=node->GetUpperNode();
+                        }
+                } //if (strcmp(field_name,"velocity")==0)
+                else if ( 
+                                (strcmp(field_name,"thickness")==0) ||
+                                (strcmp(field_name,"surface")==0)  ||
+                                (strcmp(field_name,"bed")==0)  || 
+                                (strcmp(field_name,"slopex")==0)  ||
+                                (strcmp(field_name,"slopey")==0) 
+                                ){
+
+                        /* node data: */
+                        int          doflist;
+                        int          nodedofs;
+                        double       fieldnode;
+                                
+                        GetDofList(&doflist,&numberofdofspernode);
+
+                        //initilaize node
+                        node=this;
+
+                        /*get field for this bed node: */
+                        fieldnode=field_serial[doflist];
+
+                        //go through all nodes which sit on top of this node, until we reach the surface, 
+                        //and pltg  fieldnode in field:
+                        for(;;){
+
+                                node->GetDofList(&nodedofs,&numberofdofspernode);
+                                VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES);
+
+                                if (node->IsOnSurface())break;
+                                /*get next node: */
+                                node=node->GetUpperNode();
+                        }
+                }
+                else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!"));
+
+        } //if (extrude)
+}
+                
+void  Node::UpdateNodePosition(double* thickness,double* bed){
+
+        int dof;
+
+        dof=this->GetDofList1();
+
+        if(onbed){
+
+                /*this node is on the bed, set its z position to the new bed: */
+                this->x[2]=bed[dof];
+        }
+        else if (onsurface){
+
+                /*this node is on the surface, set its z position to the new bed+ new thickness: */
+                this->x[2]=bed[dof]+thickness[dof];
+        }
+        else{
+                /*do nothing, we only update grids belonging to the exterior layers: */
+        }
+}