Changeset 10126

Timestamp:

10/06/11 09:55:17 (13 years ago)

Author:

schlegel

Message:

update for spc velocity forcing

File:

• issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (modified) (7 diffs)

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

@@ -21 +21 @@
         double  stokesreconditioning;
         bool    isstokes,ismacayealpattyn;
+   bool    spcpresent=false;
+        int Mx,Nx;
+        int My,Ny;
+        int Mz,Nz;
         double *spcvx           = NULL;
         double *spcvy           = NULL;
@@ -30 +34 @@
         double *surface         = NULL;
         double *z               = NULL;
-        
+        double *timesx=NULL;
+        double *timesy=NULL;
+        double *timesz=NULL;
+   double* values=NULL;
+
         /*Output*/
         Constraints* constraints = NULL;
@@ -61 +69 @@
         
         /*Constraints: fetch data: */
-        iomodel->FetchData(&spcvx,NULL,NULL,DiagnosticSpcvxEnum);
-        iomodel->FetchData(&spcvy,NULL,NULL,DiagnosticSpcvyEnum);
-        iomodel->FetchData(&spcvz,NULL,NULL,DiagnosticSpcvzEnum);
+        iomodel->FetchData(&spcvx,&Mx,&Nx,DiagnosticSpcvxEnum);
+        iomodel->FetchData(&spcvy,&My,&Ny,DiagnosticSpcvyEnum);
+        iomodel->FetchData(&spcvz,&Mz,&Nz,DiagnosticSpcvzEnum);
         iomodel->FetchData(&nodeonmacayeal,NULL,NULL,FlowequationBordermacayealEnum);
         if(dim==3)iomodel->FetchData(&nodeonpattyn,NULL,NULL,FlowequationBorderpattynEnum);
@@ -73 +81 @@
         /*Initialize counter: */
         count=0;
-        
+
+        /*figure out times: */
+        timesx=(double*)xmalloc(Nx*sizeof(double));
+        for(j=0;j<Nx;j++){
+                timesx[j]=spcvx[(Mx-1)*Nx+j];
+        }
+        /*unit conversion: */
+        UnitConversion(timesx,Nx,ExtToIuEnum,TimeEnum);
+        /*figure out times: */
+        timesy=(double*)xmalloc(Ny*sizeof(double));
+        for(j=0;j<Ny;j++){
+                timesy[j]=spcvy[(My-1)*Ny+j];
+        }
+        /*unit conversion: */
+        UnitConversion(timesy,Ny,ExtToIuEnum,TimeEnum);
+        /*figure out times: */
+        timesz=(double*)xmalloc(Nz*sizeof(double));
+        for(j=0;j<Nz;j++){
+                timesz[j]=spcvz[(Mz-1)*Nz+j];
+        }
+        /*unit conversion: */
+        UnitConversion(timesz,Nz,ExtToIuEnum,TimeEnum);
+
+
         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
         for (i=0;i<numberofvertices;i++){
@@ -195 +226 @@
                         /*Now add the regular spcs*/
                         else{
-                                if (!isnan(spcvx[i])){
+                                if (Mx==numberofvertices && !isnan(spcvx[i])){
                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
                                         count++;
                                 }
+                                else if (Mx==numberofvertices+1) {
+                                        /*figure out times and values: */
+                                        values=(double*)xmalloc(Nx*sizeof(double));
+                                        spcpresent=false;
+                                        for(j=0;j<Nx;j++){
+                                                values[j]=spcvx[i*Nx+j]/yts;
+                                                if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+                                        }
+
+                                        if(spcpresent){
+                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Nx,timesx,values,DiagnosticHorizAnalysisEnum));
+                                                count++;
+                                        }
+                                        xfree((void**)&values);
+                                }
                                 else if (vertices_type[i]==HutterApproximationEnum){
-                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,DiagnosticHorizAnalysisEnum)); 
-                                        count++;
-                                }
-                                
-                                if (!isnan(spcvy[i])){
+                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,DiagnosticHorizAnalysisEnum));
+                                        count++;
+                                }
+
+                                if (My==numberofvertices && !isnan(spcvy[i])){
                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
                                         count++;
                                 }
+                                else if (My==numberofvertices+1){
+                                        /*figure out times and values: */
+                                        values=(double*)xmalloc(Ny*sizeof(double));
+                                        spcpresent=false;
+                                        for(j=0;j<Ny;j++){
+                                                values[j]=spcvy[i*Ny+j]/yts;
+                                                if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+                                        }
+                                        if(spcpresent){
+                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Ny,timesy,values,DiagnosticHorizAnalysisEnum));
+                                                count++;
+                                        }
+                                        xfree((void**)&values);
+                                }
                                 else if (vertices_type[i]==HutterApproximationEnum){
-                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,DiagnosticHorizAnalysisEnum)); 
-                                        count++;
-                                }
-
-                                if (!isnan(spcvz[i]) && ((int)vertices_type[i]==StokesApproximationEnum ||  ((int)vertices_type[i]==NoneApproximationEnum))){
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
-                                        count++;
+                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,DiagnosticHorizAnalysisEnum));
+                                        count++;
+                                }
+
+                                if ((int)vertices_type[i]==StokesApproximationEnum ||  ((int)vertices_type[i]==NoneApproximationEnum)){
+                                        if (Mz==numberofvertices && !isnan(spcvz[i])){
+                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
+                                                count++;
+                                        }
+                                        else if (Mz==numberofvertices+1){
+                                                /*figure out times and values: */
+                                                values=(double*)xmalloc(Nz*sizeof(double));
+                                                spcpresent=false;
+                                                for(j=0;j<Nz;j++){
+                                                        values[j]=spcvz[i*Nz+j]/yts;
+                                                        if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+                                                }
+                                                if(spcpresent){
+                                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,Nz,timesz,values,DiagnosticHorizAnalysisEnum));
+                                                        count++;
+                                                }
+                                                xfree((void**)&values);
+                                        }
                                 }
                                 if ((int)vertices_type[i]==NoneApproximationEnum){
@@ -221 +297 @@
                                         count++;
                                 }
+
                         }
                 }
@@ -236 +313 @@
         xfree((void**)&z);
 
+        /*Free resources:*/
+        xfree((void**)&timesx);
+        xfree((void**)&timesy);
+        xfree((void**)&timesz);
+        xfree((void**)&values);
+
         /*Assign output pointer: */
         *pconstraints=constraints;