Changeset 3479 for issm/trunk

Timestamp:

04/08/10 16:49:03 (15 years ago)

Author:

Mathieu Morlighem

Message:

bug fix

Location:

issm/trunk/src

Files:

• c/objects/ElementProperties.cpp (modified) (3 diffs)
• c/objects/ElementProperties.h (modified) (2 diffs)
• c/objects/Matice.cpp (modified) (2 diffs)
• c/objects/Node.cpp (modified) (1 diff)
• c/objects/Tria.cpp (modified) (2 diffs)
• c/objects/Vertex.cpp (modified) (2 diffs)
• m/solutions/jpl/control_core.m (modified) (1 diff)
• m/solutions/jpl/diagnostic_core_linear.m (modified) (1 diff)
• m/solutions/jpl/diagnostic_core_nonlinear.m (modified) (2 diffs)
• m/solutions/jpl/thermal_core_nonlinear.m (modified) (3 diffs)
• mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (modified) (1 diff)
• mex/SystemMatrices/SystemMatrices.cpp (modified) (1 diff)
• mex/UpdateFromInputs/UpdateFromInputs.cpp (modified) (1 diff)

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

@@ -34 +34 @@
         this->geothermalflux=NULL;
         
-        this->friction_type=0;
-        this->p=0;
-        this->q=0;
-        this->shelf=0;
-        this->onbed=0;
-        this->onwater=0;
-        this->onsurface=0;
-        this->collapse=0;
-        this->thermal_steadystate=0;
+        this->friction_type=UNDEF;
+        this->p=UNDEF;
+        this->q=UNDEF;
+        this->shelf=UNDEF;
+        this->onbed=UNDEF;
+        this->onwater=UNDEF;
+        this->onsurface=UNDEF;
+        this->collapse=UNDEF;
+        this->thermal_steadystate=UNDEF;
         return;
 }
@@ -145 +145 @@
         this->numnodes=elementproperties_numnodes;
 
-        if(elementproperties_h)this->h=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->h=NULL;
-        if(elementproperties_s)this->s=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->s=NULL;
-        if(elementproperties_b)this->b=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->b=NULL;
-        if(elementproperties_k)this->k=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->k=NULL;
-        if(elementproperties_melting)this->melting=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->melting=NULL;
-        if(elementproperties_accumulation)this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->accumulation=NULL;
-        if(elementproperties_geothermalflux)this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double));
-        else this->geothermalflux=NULL;
+        this->h=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->s=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->b=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->k=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->melting=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double));
+        this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double));
 
         for(i=0;i<this->numnodes;i++){
-                if(elementproperties_h)this->h[i]=elementproperties_h[i];
-                if(elementproperties_s)this->s[i]=elementproperties_s[i];
-                if(elementproperties_b)this->b[i]=elementproperties_b[i];
-                if(elementproperties_k)this->k[i]=elementproperties_k[i];
-                if(elementproperties_melting)this->melting[i]=elementproperties_melting[i];
-                if(elementproperties_accumulation)this->accumulation[i]=elementproperties_accumulation[i];
+                if(elementproperties_h)             this->h[i]=elementproperties_h[i];
+                else                                this->h[i]=UNDEF;
+                if(elementproperties_s)             this->s[i]=elementproperties_s[i];
+                else                                this->s[i]=UNDEF;
+                if(elementproperties_b)             this->b[i]=elementproperties_b[i];
+                else                                this->b[i]=UNDEF;
+                if(elementproperties_k)             this->k[i]=elementproperties_k[i];
+                else                                this->k[i]=UNDEF;
+                if(elementproperties_melting)       this->melting[i]=elementproperties_melting[i];
+                else                                this->melting[i]=UNDEF;
+                if(elementproperties_accumulation)  this->accumulation[i]=elementproperties_accumulation[i];
+                else                                this->accumulation[i]=UNDEF;
                 if(elementproperties_geothermalflux)this->geothermalflux[i]=elementproperties_geothermalflux[i];
+                else                                this->geothermalflux[i]=UNDEF;
         }
         
@@ -179 +179 @@
         this->collapse=elementproperties_collapse;
         this->thermal_steadystate=elementproperties_thermal_steadystate;
+
+        return;
+}
+/*}}}*/
+/*FUNCTION ElementProperties Init propreties from iomodel{{{1*/
+void ElementProperties::Init(int numberofnodes, int* vertex_ids, int element_id, IoModel* iomodel){ //i is the node index
+
+        int i,k;
+
+        /*Initialize number of nodes*/
+        this->numnodes=numberofnodes;
+
+        /*Allocate node properties*/
+        if (iomodel->thickness)     this->h=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->surface)       this->s=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->bed)           this->b=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->drag)          this->k=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->melting)       this->melting=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->accumulation)  this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double));
+        if (iomodel->geothermalflux)this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double));
+
+        /*Get node properties*/
+        for(i=0;i<this->numnodes;i++){
+
+                //get index of the vertex in C indexing
+                k=vertex_ids[i]-1;
+
+                //Get properties from iomodel
+                if(iomodel->thickness)     this->h[i]=iomodel->thickness[k];
+                if(iomodel->surface)       this->s[i]=iomodel->surface[k];
+                if(iomodel->bed)           this->b[i]=iomodel->bed[k];
+                if(iomodel->drag)          this->k[i]=iomodel->drag[k];
+                if(iomodel->melting)       this->melting[i]=iomodel->melting[k];
+                if(iomodel->accumulation)  this->accumulation[i]=iomodel->accumulation[k];
+                if(iomodel->geothermalflux)this->geothermalflux[i]=iomodel->geothermalflux[k];
+        }
+
+        /*Get element properties from iomodel*/
+        if (iomodel->p)                this->p=iomodel->p[element_id-1];
+        if (iomodel->q)                this->q=iomodel->q[element_id-1];
+        if (iomodel->elementoniceshelf)this->shelf=(int)iomodel->elementoniceshelf[element_id-1];
+        if (iomodel->elementonbed)     this->onbed=(int)iomodel->elementonbed[element_id-1];
+        if (iomodel->elementonwater)   this->onwater=(bool)iomodel->elementonwater[element_id-1];
+        if (iomodel->elementonsurface) this->onsurface=(int)iomodel->elementonsurface[element_id-1];
+
+        /*Get other properties from iomodel*/
+        this->friction_type=iomodel->drag_type;
+
+        /*WARNING: collapse and thermal_steadystate are left undefined*/
 
         return;

issm/trunk/src/c/objects/ElementProperties.h

@@ -5 +5 @@
 #ifndef _ELEMENTPROPERTIES_H_
 #define  _ELEMENTPROPERTIES_H_
+
+/*indefinitions: */
+struct IoModel;
 
 class ElementProperties{
@@ -33 +36 @@
                 ElementProperties();
                 ElementProperties(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate);
-                void Init(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate);
                 ElementProperties(ElementProperties* properties);
                 ~ElementProperties();
+
+                void Init(int numnodes, int* vertex_ids, int element_id, IoModel* iomodel);
+                void Init(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate);
                 
                 void  Marshall(char** pmarshalled_dataset);

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

@@ -12 +12 @@
 #include "./Matice.h"
 #include <string.h>
+#include "../EnumDefinitions/EnumDefinitions.h"
 #include "../shared/shared.h"
 #include "../include/macros.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-
+#include "../include/typedefs.h"
                 
 /*Object constructors and destructor*/
@@ -51 +51 @@
         matice_mid=i+1;  //same as element it is created for
  
-        /*B and n: */
-        B_avg=0;
-        for(j=0;j<num_vertices;j++){
-                B_avg+=*(iomodel->B+((int)*(iomodel->elements+num_vertices*i+j)-1));
-        }
-        B_avg=B_avg/num_vertices;
-        
-        matice_B=B_avg;
-        matice_n=(double)*(iomodel->n+i);
+        /*Compute B on the element if provided*/
+        if (iomodel->B){
+                B_avg=0;
+                for(j=0;j<num_vertices;j++){
+                        B_avg+=*(iomodel->B+((int)*(iomodel->elements+num_vertices*i+j)-1));
+                }
+                B_avg=B_avg/num_vertices;
+        }
+        
+        if (iomodel->B) matice_B=B_avg;
+        else            matice_B=UNDEF;
+        if (iomodel->n) matice_n=(double)*(iomodel->n+i);
+        else            matice_n=UNDEF;
 
         this->Init(matice_mid,matice_B,matice_n);

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

@@ -803 +803 @@
 /*FUNCTION Node::UpdateFromInputs {{{2*/
 void  Node::UpdateFromInputs(void* vinputs){
-
-        ISSMERROR("not used yet!");
-        
+        
+        /*nothing updated for now*/
+
 }
 /*}}}*/

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

@@ -67 +67 @@
 
         int j;
-        int offset;
         int tria_node_ids[3];
         int tria_matice_id;
         int tria_matpar_id;
         int tria_numpar_id;
-        double tria_h[3];
-        double tria_s[3];
-        double tria_b[3];
-        double tria_k[3];
-        double tria_melting[3];
-        double tria_accumulation[3];
-        int    tria_friction_type;
-        double tria_p,tria_q;
-        int tria_shelf;
-        int tria_onwater;
-        
 
         /*id: */
@@ -99 +87 @@
         this->hmatpar.Init(&tria_matpar_id,1);
         this->hnumpar.Init(&tria_numpar_id,1);
-
-        /*properties: */
-        for(j=0;j<3;j++){
-
-                offset=(int)*(iomodel->elements+3*i+j)-1; //get index of this vertex, in C numbering. 
-                tria_h[j]= *(iomodel->thickness+offset); 
-                tria_s[j]= *(iomodel->surface+offset); 
-                tria_b[j]=*(iomodel->bed+offset);
-                tria_k[j]=*(iomodel->drag+offset);
-                tria_melting[j]=*(iomodel->melting+offset);
-                tria_accumulation[j]=*(iomodel->accumulation+offset);
-        }
-        tria_friction_type=(int)iomodel->drag_type;
-        tria_p=iomodel->p[i];
-        tria_q=iomodel->q[i];
-        tria_shelf=(int)*(iomodel->elementoniceshelf+i);
-        tria_onwater=(bool)*(iomodel->elementonwater+i);
-        
-        this->properties.Init(3,tria_h, tria_s, tria_b, tria_k, tria_melting, tria_accumulation, NULL,
-                        tria_friction_type, tria_p, tria_q, tria_shelf, UNDEF,tria_onwater, UNDEF,UNDEF,UNDEF);
-
+        this->properties.Init(3,tria_node_ids,this->id,iomodel);
 
 }

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

@@ -188 +188 @@
         Vertex* vertex=NULL;
         int dof[1]={0};
-        double coord[3];
 
         vertex=this;
-
-        coord[0]=this->x; coord[1]=this->y; coord[2]=this->z;
 
         /*Recover parameter inputs: */
@@ -198 +195 @@
 
         /*Update internal data if inputs holds new values: */
-        inputs->Recover("x",&coord[0],1,dof,1,(void**)&vertex);
-        inputs->Recover("y",&coord[1],1,dof,1,(void**)&vertex);
-        inputs->Recover("z",&coord[2],1,dof,1,(void**)&vertex);
-        
-        ISSMERROR("not supported yet!");
+        inputs->Recover("x",&x,1,dof,1,(void**)&vertex);
+        inputs->Recover("y",&y,1,dof,1,(void**)&vertex);
+        inputs->Recover("z",&z,1,dof,1,(void**)&vertex);
         
 }

issm/trunk/src/m/solutions/jpl/control_core.m

@@ -45 +45 @@
 
         %Update inputs in datasets
-        [model.elements,model.nodes,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.loads,model.materials,model.parameters,inputs);
+        [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters,inputs);
 
         displaystring(verbose,'\n%s',['      computing gradJ...']);

issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m

@@ -9 +9 @@
 
         %Update inputs in datasets
-        [m.elements,m.nodes, m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, m.loads,m.materials,m.parameters,inputs);
+        [m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs);
         
         %system matrices
-        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-        [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+        [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
         
         %Reduce tangent matrix from g size to f size

issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m

@@ -31 +31 @@
                 
                 %Update inputs in datasets
-                [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs);
+                [m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters,inputs);
                 
                 %system matrices 
-                [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
         
                 %penalties
-                [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
                 %Reduce tangent matrix from g size to f size
@@ -84 +84 @@
                 inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
                 m.parameters.kflag=1; m.parameters.pflag=0; 
-                [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
                 [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
                 varargout(1)={K_ff};

issm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m

@@ -23 +23 @@
 
                 %Update inputs in datasets
-                [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs);
+                [m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters,inputs);
 
                 %system matrices 
@@ -29 +29 @@
                         if count==1
                                 displaystring(m.parameters.verbose,'%s',['   system matrices']);
-                                [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                                [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
                         end
                         displaystring(m.parameters.verbose,'%s',['   penalty system matrices']);
-                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.verticesloads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
                 else
                         displaystring(m.parameters.verbose,'%s',['   system matrices']);
-                        [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                        [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
                         displaystring(m.parameters.verbose,'%s',['   penalty system matrices']);
-                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                        [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
                 end
 
@@ -57 +57 @@
                 inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
                 displaystring(m.parameters.verbose,'%s',['   update inputs']);
-                [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs);
+                [m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters,inputs);
         
                 %penalty constraints

issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp

@@ -79 +79 @@
 {
         _printf_("\n");
-        _printf_("   usage: [Kgg,pg] = %s(Kggin,pgin,eleemnts,nodes,loads,materials,params,inputs,analysis_type,sub_analysis_type);\n",__FUNCT__);
+        _printf_("   usage: [Kgg,pg] = %s(Kggin,pgin,elements,nodes,vertices,loads,materials,params,inputs,analysis_type,sub_analysis_type);\n",__FUNCT__);
         _printf_("\n");
 }

issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

@@ -81 +81 @@
 {
         _printf_("\n");
-        _printf_("   usage: [Kgg,pg] = %s(eleemnts,nodes,loads,materials,params,inputs,analysis_type);\n",__FUNCT__);
+        _printf_("   usage: [Kgg,pg] = %s(elements,nodes,vertices,loads,materials,params,inputs,analysis_type);\n",__FUNCT__);
         _printf_("\n");
 }

issm/trunk/src/mex/UpdateFromInputs/UpdateFromInputs.cpp

@@ -64 +64 @@
 {
         _printf_("\n");
-        _printf_("   usage: [elements,nodes,loads, materials, parameters] = %s(elements,nodes,loads, materials,parameters,inputs);\n",__FUNCT__);
+        _printf_("   usage: [elements,nodes,vertices,loads,materials,parameters] = %s(elements,nodes,vertices,loads,materials,parameters,inputs);\n",__FUNCT__);
         _printf_("\n");
 }