Changeset 17537

Timestamp:

03/25/14 11:39:00 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: creating Dstar matrix

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -3419 +3419 @@
         /*Initialize vectors*/
         ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
-        IssmDouble*    Dstar   = xNew<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
+        IssmDouble*    Dstar   = xNewZeroInit<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
         IssmDouble*    tau     = xNew<IssmDouble>(tausize*tnumnodes);
         IssmDouble*    d       = xNew<IssmDouble>(tausize*tnumnodes);
@@ -3478 +3478 @@
 
                 /*Create Dstar*/
+                /*In dim = 2
+                 *
+                 *       <----------------- tausize ---------------> x tnumnodes
+                 *       |  gamma_ij^x         0       gamma_ij^y  | ^
+                 * Dij = |                                         | dim
+                 *       |     0          gamma_ij^y   gamma_ij^x  | v
+                 *                                                   x
+                 *                                                   vnumnodes
+                 *
+                 *In dim = 3
+                 *
+                 *       |  gamma_ij^x         0          0         gamma_ij^y  gamma_ij^z      0      |
+                 *       |                                                                             |
+                 * Dij = |     0          gamma_ij^y      0         gamma_ij^x     0        gamma_ij^z |
+                 *       |                                                                             |
+                 *       |     0               0      gamma_ij^z        0       gamma_ij^x  gamma_ij^y |
+                 *
+                 * gamma_ij^x = zeta_j dphi_i/dx
+                 *
+                 * where:
+                 *   - zeta_j is the nodal function for the j^th node of the tensor (P1DG)
+                 *   - phi_i  is the nodal function for the i^th node of the velocity (P2)*/
                 element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
                 element->NodalFunctionsTensor(tbasis,gauss);
-
-                _error_("STOP");
+                if(dim==2){
+                        for(int i=0;i<vnumnodes;i++){
+                                for(int j=0;j<tnumnodes;j++){
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] = tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[1*vnumnodes+i];
+
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] = tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[0*vnumnodes+i];
+                                }
+                        }
+                }
+                else{
+                        for(int i=0;i<vnumnodes;i++){
+                                for(int j=0;j<tnumnodes;j++){
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] = tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] = tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] = tbasis[j]*vdbasis[2*vnumnodes+i];
+
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] = tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] = tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] = tbasis[j]*vdbasis[2*vnumnodes+i];
+
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[2*vnumnodes+i];
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] = tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] = tbasis[j]*vdbasis[1*vnumnodes+i];
+                                }
+                        }
+                }
         }