Changeset 17628

Timestamp:

04/02/14 16:11:41 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing Dstar

File:

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

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

@@ -3503 +3503 @@
                         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];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
                                 }
                         }
@@ -3514 +3514 @@
                         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];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+                                        Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
                                 }
                         }
                 }
-
-                /*contribution -Dstar tau*/
-                for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = -gauss->weight*Jdet;
-                TripleMultiply(Dstar,dim*tnumnodes,tausize*tnumnodes,0,
-                                        D,tausize*tnumnodes,tausize*tnumnodes,0,
-                                        tau,tausize*tnumnodes,1,0,
-                                        &pe->values[0],1);
-
-                /*contribution + r Dstar d*/
-                for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = +r*gauss->weight*Jdet;
-                TripleMultiply(Dstar,dim*tnumnodes,tausize*tnumnodes,0,
-                                        D,tausize*tnumnodes,tausize*tnumnodes,0,
-                                        d,tausize*tnumnodes,1,0,
-                                        &pe->values[0],1);
-        }
+        }
+
+        /*contribution -Dstar tau*/
+        for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = -1.;
+        TripleMultiply(Dstar,dim*vnumnodes,tausize*tnumnodes,0,
+                                D,tausize*tnumnodes,tausize*tnumnodes,0,
+                                tau,tausize*tnumnodes,1,0,
+                                &pe->values[0],1);
+
+        /*contribution + r Dstar d*/
+        for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = +r;
+        TripleMultiply(Dstar,dim*vnumnodes,tausize*tnumnodes,0,
+                                D,tausize*tnumnodes,tausize*tnumnodes,0,
+                                d,tausize*tnumnodes,1,0,
+                                &pe->values[0],1);
 
         /*Transform coordinate system*/