source: issm/oecreview/Archive/24684-25833/ISSM-25262-25263.diff@ 25834

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

CHG: added 24684-25833
File size: 10.7 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
2===================================================================
3--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 25262)
4+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 25263)
5@@ -29,8 +29,6 @@
6 ElementMatrix* CreateJacobianMatrix(Element* element);
7 ElementMatrix* CreateKMatrix(Element* element);
8 ElementVector* CreatePVector(Element* element);
9- void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
10- void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
11 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
12 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
13 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
14Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
15===================================================================
16--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 25262)
17+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 25263)
18@@ -453,14 +453,12 @@
19 /*Initialize Element vector*/
20 ElementMatrix* Ke = element->NewElementMatrix();
21 IssmDouble* basis = xNew<IssmDouble>(numnodes);
22- IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
23- IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
24+ IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
25 IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
26
27 /*Retrieve all inputs and parameters*/
28 element->GetVerticesCoordinates(&xyz_list);
29 element->FindParam(&dt,TimesteppingTimeStepEnum);
30- //printf("dt %f\n", dt);
31 element->FindParam(&stabilization,DamageStabilizationEnum);
32 Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
33 Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
34@@ -478,6 +476,7 @@
35
36 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
37 element->NodalFunctions(basis,gauss);
38+ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
39
40 vx_input->GetInputValue(&vx,gauss);
41 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
42@@ -489,38 +488,36 @@
43 vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
44 }
45
46+ /*Transient term*/
47 D_scalar=gauss->weight*Jdet;
48- TripleMultiply(basis,1,numnodes,1,
49- &D_scalar,1,1,0,
50- basis,1,numnodes,0,
51- &Ke->values[0],1);
52+ for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
53
54- GetB(B,element,dim,xyz_list,gauss);
55- GetBprime(Bprime,element,dim,xyz_list,gauss);
56-
57 dvxdx=dvx[0];
58 dvydy=dvy[1];
59- if(dim==3) dvzdz=dvz[2];
60 D_scalar=dt*gauss->weight*Jdet;
61+ if(dim==2){
62+ for(int i=0;i<numnodes;i++){
63+ for(int j=0;j<numnodes;j++){
64+ /*\phi_i \phi_j \nabla\cdot v*/
65+ Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
66+ /*\phi_i v\cdot\nabla\phi_j*/
67+ Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
68+ }
69+ }
70+ }
71+ else{/*3D*/
72+ _assert_(dim==3);
73+ dvzdz=dvz[2];
74+ for(int i=0;i<numnodes;i++){
75+ for(int j=0;j<numnodes;j++){
76+ /*\phi_i \phi_j \nabla\cdot v*/
77+ Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy+dvzdz);
78+ /*\phi_i v\cdot\nabla\phi_j*/
79+ Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j] + vz*dbasis[2*numnodes+j]);
80+ }
81+ }
82+ }
83
84- D[0*dim+0]=D_scalar*dvxdx;
85- D[1*dim+1]=D_scalar*dvydy;
86- if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
87-
88- TripleMultiply(B,dim,numnodes,1,
89- D,dim,dim,0,
90- B,dim,numnodes,0,
91- &Ke->values[0],1);
92-
93- D[0*dim+0]=D_scalar*vx;
94- D[1*dim+1]=D_scalar*vy;
95- if(dim==3) D[2*dim+2]=D_scalar*vz;
96-
97- TripleMultiply(B,dim,numnodes,1,
98- D,dim,dim,0,
99- Bprime,dim,numnodes,0,
100- &Ke->values[0],1);
101-
102 if(stabilization==2){
103 if(dim==3){
104 vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8;
105@@ -565,6 +562,14 @@
106 D[1*dim+0]=D_scalar*D[1*dim+0];
107 D[0*dim+1]=D_scalar*D[0*dim+1];
108 D[1*dim+1]=D_scalar*D[1*dim+1];
109+ for(int i=0;i<numnodes;i++){
110+ for(int j=0;j<numnodes;j++){
111+ Ke->values[i*numnodes+j] += (
112+ dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
113+ dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
114+ );
115+ }
116+ }
117 }
118 else if(dim==3){
119 D[0*dim+0]=D_scalar*D[0*dim+0];
120@@ -576,20 +581,23 @@
121 D[0*dim+2]=D_scalar*D[0*dim+2];
122 D[1*dim+2]=D_scalar*D[1*dim+2];
123 D[2*dim+2]=D_scalar*D[2*dim+2];
124+ for(int i=0;i<numnodes;i++){
125+ for(int j=0;j<numnodes;j++){
126+ Ke->values[i*numnodes+j] += (
127+ dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j] + D[0*dim+2]*dbasis[2*numnodes+j]) +
128+ dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j] + D[1*dim+2]*dbasis[2*numnodes+j]) +
129+ dbasis[2*numnodes+i] *(D[2*dim+0]*dbasis[0*numnodes+j] + D[2*dim+1]*dbasis[1*numnodes+j] + D[2*dim+2]*dbasis[2*numnodes+j])
130+ );
131+ }
132+ }
133 }
134- TripleMultiply(Bprime,dim,numnodes,1,
135- D,dim,dim,0,
136- Bprime,dim,numnodes,0,
137- &Ke->values[0],1);
138 }
139-
140 }
141
142 /*Clean up and return*/
143 xDelete<IssmDouble>(xyz_list);
144 xDelete<IssmDouble>(basis);
145- xDelete<IssmDouble>(B);
146- xDelete<IssmDouble>(Bprime);
147+ xDelete<IssmDouble>(dbasis);
148 xDelete<IssmDouble>(D);
149 delete gauss;
150 return Ke;
151@@ -665,63 +673,6 @@
152 delete gauss;
153 return pe;
154 }/*}}}*/
155-void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
156- /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
157- * For node i, Bi can be expressed in the actual coordinate system
158- * by:
159- * Bi=[ N ]
160- * [ N ]
161- * where N is the finiteelement function for node i.
162- *
163- * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
164- */
165-
166- /*Fetch number of nodes for this finite element*/
167- int numnodes = element->GetNumberOfNodes();
168-
169- /*Get nodal functions*/
170- IssmDouble* basis=xNew<IssmDouble>(numnodes);
171- element->NodalFunctions(basis,gauss);
172-
173- /*Build B: */
174- for(int i=0;i<numnodes;i++){
175- for(int j=0;j<dim;j++){
176- B[numnodes*j+i] = basis[i];
177- }
178- }
179-
180- /*Clean-up*/
181- xDelete<IssmDouble>(basis);
182-}/*}}}*/
183-void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
184- /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
185- * For node i, Bi' can be expressed in the actual coordinate system
186- * by:
187- * Bi_prime=[ dN/dx ]
188- * [ dN/dy ]
189- * where N is the finiteelement function for node i.
190- *
191- * We assume B' has been allocated already, of size: 3x(2*numnodes)
192- */
193-
194- /*Fetch number of nodes for this finite element*/
195- int numnodes = element->GetNumberOfNodes();
196-
197- /*Get nodal functions derivatives*/
198- IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
199- element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
200-
201- /*Build B': */
202- for(int i=0;i<numnodes;i++){
203- for(int j=0;j<dim;j++){
204- Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
205- }
206- }
207-
208- /*Clean-up*/
209- xDelete<IssmDouble>(dbasis);
210-
211-}/*}}}*/
212 void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
213 element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
214 }/*}}}*/
215@@ -787,9 +738,8 @@
216
217 /*Initialize Element vector and other vectors*/
218 ElementMatrix* Ke = element->NewElementMatrix();
219- IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
220- IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
221- IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
222+ IssmDouble* basis = xNew<IssmDouble>(numnodes);
223+ IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
224
225 /*Retrieve all inputs and parameters*/
226 element->GetVerticesCoordinates(&xyz_list);
227@@ -802,26 +752,24 @@
228 gauss->GaussPoint(ig);
229
230 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
231- GetB(B,element,dim,xyz_list,gauss);
232- GetBprime(Bprime,element,dim,xyz_list,gauss);
233+ element->NodalFunctions(basis,gauss);
234+ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
235+
236 vxaverage_input->GetInputValue(&vx,gauss);
237 vyaverage_input->GetInputValue(&vy,gauss);
238
239- D[0*dim+0] = -gauss->weight*vx*Jdet;
240- D[1*dim+1] = -gauss->weight*vy*Jdet;
241-
242- TripleMultiply(B,dim,numnodes,1,
243- D,dim,dim,0,
244- Bprime,dim,numnodes,0,
245- &Ke->values[0],1);
246-
247+ for(int i=0;i<numnodes;i++){
248+ for(int j=0;j<numnodes;j++){
249+ /*\phi_i v\cdot\nabla\phi_j*/
250+ Ke->values[i*numnodes+j] += -gauss->weight*Jdet*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
251+ }
252+ }
253 }
254
255 /*Clean up and return*/
256 xDelete<IssmDouble>(xyz_list);
257- xDelete<IssmDouble>(B);
258- xDelete<IssmDouble>(Bprime);
259- xDelete<IssmDouble>(D);
260+ xDelete<IssmDouble>(basis);
261+ xDelete<IssmDouble>(dbasis);
262 delete gauss;
263 return Ke;
264 }/*}}}*/
265@@ -853,10 +801,7 @@
266 element->NodalFunctions(basis,gauss);
267
268 D=gauss->weight*Jdet;
269- TripleMultiply(basis,1,numnodes,1,
270- &D,1,1,0,
271- basis,1,numnodes,0,
272- &Me->values[0],1);
273+ for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Me->values[i*numnodes+j] += D*basis[i]*basis[j];
274 }
275
276 /*Clean up and return*/
277Index: ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp
278===================================================================
279--- ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25262)
280+++ ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25263)
281@@ -107,11 +107,7 @@
282 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
283 basalelement->NodalFunctions(basis,gauss);
284 D=gauss->weight*Jdet;
285-
286- TripleMultiply(basis,1,numnodes,1,
287- &D,1,1,0,
288- basis,1,numnodes,0,
289- &Ke->values[0],1);
290+ for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
291 }
292
293 /*Clean up and return*/
294Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
295===================================================================
296--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25262)
297+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25263)
298@@ -336,11 +336,7 @@
299
300 /*Transient term*/
301 D_scalar=gauss->weight*Jdet;
302- for(int i=0;i<numnodes;i++){
303- for(int j=0;j<numnodes;j++){
304- Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
305- }
306- }
307+ for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
308
309 /*Advection terms*/
310 vxaverage_input->GetInputValue(&vx,gauss);
Note: See TracBrowser for help on using the repository browser.