Changeset 3556

Timestamp:

04/16/10 11:34:43 (15 years ago)

Author:

Mathieu Morlighem

Message:

fixed prognostic2

Location:

issm/trunk/src

Files:

• c/DataSet/DataSet.cpp (modified) (2 diffs)
• c/DataSet/DataSet.h (modified) (1 diff)
• c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.cpp (modified) (2 diffs)
• c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.h (modified) (1 diff)
• c/objects/Node.cpp (modified) (1 diff)
• c/objects/Node.h (modified) (1 diff)
• c/objects/Tria.cpp (modified) (1 diff)
• c/parallel/ProcessResults.cpp (modified) (2 diffs)
• c/parallel/prognostic2_core.cpp (modified) (1 diff)
• c/shared/Elements/elements.h (modified) (1 diff)
• m/solutions/jpl/prognostic2_core.m (modified) (1 diff)
• mex/FieldAverageOntoVertices/FieldAverageOntoVertices.cpp (modified) (4 diffs)
• mex/FieldAverageOntoVertices/FieldAverageOntoVertices.h (modified) (2 diffs)

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

@@ -1073 +1073 @@
 /*}}}*/
 /*FUNCTION DataSet::FieldAverageOntoVertices{{{1*/
-void  DataSet::FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname){
+void DataSet::FieldAverageOntoVertices(Vec fieldsum,Vec connectivity,double* field){
 
         vector<Object*>::iterator object;
@@ -1082 +1082 @@
                 if((*object)->Enum()==NodeEnum()){
                         node=(Node*)(*object);
-                        node->FieldAverageOntoVertices(field,field_serial,fieldname);
+                        node->FieldAverageOntoVertices(fieldsum,connectivity,field);
                 }
 

issm/trunk/src/c/DataSet/DataSet.h

@@ -88 +88 @@
                 void  CostFunction(double* pJ, void* inputs,int analysis_type,int sub_analysis_type);
                 void  SurfaceArea(double* pS, void* inputs,int analysis_type,int sub_analysis_type);
-                void  FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname);
+                void  FieldAverageOntoVertices(Vec fieldsum,Vec connectivity,double* field);
                 void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
                 int   DeleteObject(Object* object);

issm/trunk/src/c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.cpp

@@ -10 +10 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-void FieldAverageOntoVerticesx( Vec field, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters,char* fieldname){
+void FieldAverageOntoVerticesx(Vec* pfield, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters){
 
+        extern int my_rank;
+
+        /*Intermediary*/
+        int     numberofvertices,found;
+        Vec     connectivity=NULL;
+        Vec     fieldsum=NULL;
         double* field_serial=NULL;
+
+        /*Output*/
+        Vec field=NULL;
+
+        /*Initioalize intermediary*/
+        found=parameters->FindParam(&numberofvertices,"numberofvertices");
+        connectivity=NewVec(numberofvertices);
+        fieldsum    =NewVec(numberofvertices);
 
         /*First, get elements and nodes configured: */
@@ -20 +34 @@
 
         /*Serialize field: */
-        VecToMPISerial(&field_serial,field);
+        VecToMPISerial(&field_serial,*pfield);
 
-        /*Depth average field, onto base of mesh: */
-        nodes->FieldAverageOntoVertices(field,field_serial,fieldname);
+        /*Compute summ of connectivity and summ of values: */
+        nodes->FieldAverageOntoVertices(fieldsum,connectivity,field_serial);
 
-        /*Assemble vector: */
-        VecAssemblyBegin(field);
-        VecAssemblyEnd(field);
+        /*Assemble vectors: */
+        VecAssemblyBegin(connectivity);
+        VecAssemblyEnd(connectivity);
+        VecAssemblyBegin(fieldsum);
+        VecAssemblyEnd(fieldsum);
+
+        /*Initialize output*/
+        field=NewVec(numberofvertices);
+
+        /*Create Field = Fieldsum./connectivity */
+        VecPointwiseDivide(field,fieldsum,connectivity);
+
+        /*Assign output pointer*/
+        VecFree(pfield);
+        *pfield=field;
 
         /*Free ressources:*/
         xfree((void**)&field_serial);
+        VecFree(&connectivity);
+        VecFree(&fieldsum);
 
 }

issm/trunk/src/c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.h

@@ -9 +9 @@
 
 /* local prototypes: */
-void FieldAverageOntoVerticesx( Vec field, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters,char* fieldname);
+void FieldAverageOntoVerticesx(Vec* field, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters);
 
 #endif  /* _FIELDDEPTHAVERAGEX_H */

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

@@ -484 +484 @@
 /*}}}*/
 /*FUNCTION Node::FieldAverageOntoVertices{{{2*/
-void  Node::FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname){
-
-        ISSMERROR("NOT SUPPORTED YET"); // we need the vertex number !
+void  Node::FieldAverageOntoVertices(Vec fieldsum,Vec connectivity,double* field){
+
+        /*Intermediary*/
+        int vertexdof;
+        int index;
+        double value;
+
+        /*Skip if clone*/
+        if (indexing.clone==1) return;
+
+        /*Check dofpernode*/
+        if (indexing.numberofdofs!=1) ISSMERROR("only one dof suppoerted yet");
+
+        /*Get value at node*/
+        index=indexing.doflist[0];
+        value=field[index];
+        printf("value[%i]=%g\n",id,value);
+
+        /*Add values to the two vectors*/
+        vertexdof=this->GetVertexDof();
+        VecSetValue(fieldsum,    vertexdof,value,ADD_VALUES);
+        VecSetValue(connectivity,vertexdof,1.0,  ADD_VALUES);
 }
 /*}}}*/

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

@@ -83 +83 @@
                 int   IsOnSurface();
                 void  FreezeDof(int dof);
-                void  FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname);
+                void  FieldAverageOntoVertices(Vec fieldsum,Vec connectivity,double* field);
                 void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
                 int   IsOnShelf();

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

@@ -1959 +1959 @@
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-#ifdef _DEBUGELEMENTS_
-        if(my_rank==RANK && id==ELID){ 
-                printf("      Ke_gg erms:\n");
-                for( i=0; i<numdof; i++){
-                        for (j=0;j<numdof;j++){
-                                printf("%g ",Ke_gg[i][j]);
-                        }
-                        printf("\n");
-                }
-                printf("      Ke_gg row_indices:\n");
-                for( i=0; i<numdof; i++){
-                        printf("%i ",doflist[i]);
-                }
-
-        }
-#endif
 
 cleanup_and_return: 

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

@@ -106 +106 @@
         double* v_g_serial=NULL;
 
-        int numberofnodes,numberofelements;
+        int numberofnodes,numberofvertices,numberofelements;
 
         /*Initialize new results: */
@@ -350 +350 @@
                 }
                 else if(strcmp(result->GetFieldName(),"h_g")==0){
-                        /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
+                        /*easy, h_g is of size numberofvertices, on 1 dof, just repartition: */
                         result->GetField(&h_g);
                         VecToMPISerial(&h_g_serial,h_g);
-                        fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
+                        fem_p->parameters->FindParam(&numberofvertices,"numberofvertices");
                         VecToMPISerial(&partition,fem_p->partition->vector);
 
-                        thickness=(double*)xmalloc(numberofnodes*sizeof(double));
-
-                        for(i=0;i<numberofnodes;i++){
+                        thickness=(double*)xmalloc(numberofvertices*sizeof(double));
+
+                        for(i=0;i<numberofvertices;i++){
                                 thickness[i]=h_g_serial[(int)partition[i]];
                         }
                         
                         /*Ok, add pressure to newresults: */
-                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
+                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofvertices);
                         newresults->AddObject(newresult);
 

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

@@ -56 +56 @@
         diagnostic_core_linear(&h_g,fem_p,inputs,Prognostic2AnalysisEnum(),NoneAnalysisEnum());
 
-        _printf_("extrude computed thickness on all layers:\n");
-        //FieldAverageOntoVertices(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+        _printf_("Averaging over vertices:\n");
+        FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
+
+        //_printf_("extrude computed thickness on all layers:\n");
+        //FieldExtrudex(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
 
         /*Plug results into output dataset: */

issm/trunk/src/c/shared/Elements/elements.h

@@ -18 +18 @@
                 printf("   [ ");
                 for (j=0;j<cols;j++){
-                        printf(" %10g ",array[i*cols+j]);
+                        printf(" %10.10g ",array[i*cols+j]);
+                }  
+                printf(" ]\n");
+        }  
+        printf("\n");
+}
+inline void printarray(int* array,int lines,int cols){
+        int i,j;
+        printf("\n");
+        for (i=0;i<lines;i++){  
+                printf("   [ ");
+                for (j=0;j<cols;j++){
+                        printf(" %i",array[i*cols+j]);
                 }  
                 printf(" ]\n");

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

@@ -25 +25 @@
         results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
 
+        displaystring(m.parameters.verbose,'\n%s',['averaging over vertices']);
+        results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
+
         displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
         %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');

issm/trunk/src/mex/FieldAverageOntoVertices/FieldAverageOntoVertices.cpp

@@ -18 +18 @@
         DataSet* parameters=NULL;
         Vec      field=NULL;
-        char*    fieldname=NULL;
 
         /*Boot module: */
@@ -34 +33 @@
         FetchParams(&parameters,PARAMETERS);
         FetchData(&field,FIELD);
-        FetchData(&fieldname,FIELDNAME);
 
         /*!Call core code: */
-        FieldAverageOntoVerticesx(field,elements,nodes,vertices,loads,materials,parameters,fieldname);
+        FieldAverageOntoVerticesx(&field,elements,nodes,vertices,loads,materials,parameters);
 
         /*write output : */
@@ -50 +48 @@
         delete parameters;
         VecFree(&field);
-        xfree((void**)&fieldname);
         
         /*end module: */
@@ -60 +57 @@
 {
         _printf_("\n");
-        _printf_("   usage: [field] = %s(elements, nodes,loads, materials, parameters,field,fieldname);\n",__FUNCT__);
+        _printf_("   usage: [field] = %s(elements, nodes,loads, materials, parameters,field);\n",__FUNCT__);
         _printf_("\n");
 }

issm/trunk/src/mex/FieldAverageOntoVertices/FieldAverageOntoVertices.h

@@ -24 +24 @@
 #define PARAMETERS (mxArray*)prhs[5]
 #define FIELD (mxArray*)prhs[6]
-#define FIELDNAME (mxArray*)prhs[7]
 
 /* serial output macros: */
@@ -33 +32 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  8
+#define NRHS  7