source: issm/oecreview/Archive/24684-25833/ISSM-24917-24918.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 2.0 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24917)
4+++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24918)
5@@ -590,7 +590,6 @@
6 ElementMatrix* Ke = element->NewElementMatrix();
7 IssmDouble* basis = xNew<IssmDouble>(numnodes);
8 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
9- IssmDouble* Bprime = xNew<IssmDouble>(3*numnodes);
10 IssmDouble K[3][3];
11
12 /*Retrieve all inputs and parameters*/
13@@ -666,11 +665,15 @@
14 K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
15 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
16
17- GetBAdvecprime(Bprime,element,xyz_list,gauss);
18- TripleMultiply(Bprime,3,numnodes,1,
19- &K[0][0],3,3,0,
20- Bprime,3,numnodes,0,
21- &Ke->values[0],1);
22+ for(int i=0;i<numnodes;i++){
23+ for(int j=0;j<numnodes;j++){
24+ Ke->values[i*numnodes+j] += (
25+ dbasis[0*numnodes+i] *(K[0][0]*dbasis[0*numnodes+j] + K[0][1]*dbasis[1*numnodes+j]+ K[0][2]*dbasis[2*numnodes+j]) +
26+ dbasis[1*numnodes+i] *(K[1][0]*dbasis[0*numnodes+j] + K[1][1]*dbasis[1*numnodes+j]+ K[1][2]*dbasis[2*numnodes+j]) +
27+ dbasis[2*numnodes+i] *(K[2][0]*dbasis[0*numnodes+j] + K[2][1]*dbasis[1*numnodes+j]+ K[2][2]*dbasis[2*numnodes+j])
28+ );
29+ }
30+ }
31 }
32 /*SUPG*/
33 else if(stabilization==2){
34@@ -714,7 +717,6 @@
35 xDelete<IssmDouble>(xyz_list);
36 xDelete<IssmDouble>(basis);
37 xDelete<IssmDouble>(dbasis);
38- xDelete<IssmDouble>(Bprime);
39 delete gauss;
40 return Ke;
41 }/*}}}*/
42@@ -757,11 +759,7 @@
43
44 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
45 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
46- TripleMultiply(basis,numnodes,1,0,
47- &D,1,1,0,
48- basis,1,numnodes,0,
49- &Ke->values[0],1);
50-
51+ for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
52 }
53
54 /*Clean up and return*/
Note: See TracBrowser for help on using the repository browser.