1 | /*!\file PentaVertexInput.c
|
---|
2 | * \brief: implementation of the PentaVertexInput object
|
---|
3 | */
|
---|
4 |
|
---|
5 | #ifdef HAVE_CONFIG_H
|
---|
6 | #include "config.h"
|
---|
7 | #else
|
---|
8 | #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
|
---|
9 | #endif
|
---|
10 |
|
---|
11 | #include "stdio.h"
|
---|
12 | #include <string.h>
|
---|
13 | #include "../objects.h"
|
---|
14 | #include "../../EnumDefinitions/EnumDefinitions.h"
|
---|
15 | #include "../../shared/shared.h"
|
---|
16 | #include "../../Container/Container.h"
|
---|
17 | #include "../../include/include.h"
|
---|
18 |
|
---|
19 | /*PentaVertexInput constructors and destructor*/
|
---|
20 | /*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
|
---|
21 | PentaVertexInput::PentaVertexInput(){
|
---|
22 | return;
|
---|
23 | }
|
---|
24 | /*}}}*/
|
---|
25 | /*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
|
---|
26 | PentaVertexInput::PentaVertexInput(int in_enum_type,double* in_values)
|
---|
27 | :PentaRef(1)
|
---|
28 | {
|
---|
29 |
|
---|
30 | /*Set PentaRef*/
|
---|
31 | this->SetElementType(P1Enum,0);
|
---|
32 | this->element_type=P1Enum;
|
---|
33 |
|
---|
34 | enum_type=in_enum_type;
|
---|
35 | values[0]=in_values[0];
|
---|
36 | values[1]=in_values[1];
|
---|
37 | values[2]=in_values[2];
|
---|
38 | values[3]=in_values[3];
|
---|
39 | values[4]=in_values[4];
|
---|
40 | values[5]=in_values[5];
|
---|
41 | }
|
---|
42 | /*}}}*/
|
---|
43 | /*FUNCTION PentaVertexInput::~PentaVertexInput(){{{1*/
|
---|
44 | PentaVertexInput::~PentaVertexInput(){
|
---|
45 | return;
|
---|
46 | }
|
---|
47 | /*}}}*/
|
---|
48 |
|
---|
49 | /*Object virtual functions definitions:*/
|
---|
50 | /*FUNCTION PentaVertexInput::Echo {{{1*/
|
---|
51 | void PentaVertexInput::Echo(void){
|
---|
52 | this->DeepEcho();
|
---|
53 | }
|
---|
54 | /*}}}*/
|
---|
55 | /*FUNCTION PentaVertexInput::DeepEcho{{{1*/
|
---|
56 | void PentaVertexInput::DeepEcho(void){
|
---|
57 |
|
---|
58 | printf("PentaVertexInput:\n");
|
---|
59 | printf(" enum: %i (%s)\n",this->enum_type,EnumAsString(this->enum_type));
|
---|
60 | printf(" values: [%g %g %g %g %g %g]\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]);
|
---|
61 | }
|
---|
62 | /*}}}*/
|
---|
63 | /*FUNCTION PentaVertexInput::Id{{{1*/
|
---|
64 | int PentaVertexInput::Id(void){ return -1; }
|
---|
65 | /*}}}*/
|
---|
66 | /*FUNCTION PentaVertexInput::MyRank{{{1*/
|
---|
67 | int PentaVertexInput::MyRank(void){
|
---|
68 | extern int my_rank;
|
---|
69 | return my_rank;
|
---|
70 | }
|
---|
71 | /*}}}*/
|
---|
72 | /*FUNCTION PentaVertexInput::Marshall{{{1*/
|
---|
73 | void PentaVertexInput::Marshall(char** pmarshalled_dataset){
|
---|
74 |
|
---|
75 | char* marshalled_dataset=NULL;
|
---|
76 | int enum_value=0;
|
---|
77 |
|
---|
78 | /*recover marshalled_dataset: */
|
---|
79 | marshalled_dataset=*pmarshalled_dataset;
|
---|
80 |
|
---|
81 | /*get enum value of PentaVertexInput: */
|
---|
82 | enum_value=PentaVertexInputEnum;
|
---|
83 |
|
---|
84 | /*marshall enum: */
|
---|
85 | memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
|
---|
86 |
|
---|
87 | /*marshall PentaVertexInput data: */
|
---|
88 | memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
|
---|
89 | memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);
|
---|
90 |
|
---|
91 | *pmarshalled_dataset=marshalled_dataset;
|
---|
92 | }
|
---|
93 | /*}}}*/
|
---|
94 | /*FUNCTION PentaVertexInput::MarshallSize{{{1*/
|
---|
95 | int PentaVertexInput::MarshallSize(){
|
---|
96 |
|
---|
97 | return sizeof(values)+
|
---|
98 | +sizeof(enum_type)+
|
---|
99 | +sizeof(int); //sizeof(int) for enum value
|
---|
100 | }
|
---|
101 | /*}}}*/
|
---|
102 | /*FUNCTION PentaVertexInput::Demarshall{{{1*/
|
---|
103 | void PentaVertexInput::Demarshall(char** pmarshalled_dataset){
|
---|
104 |
|
---|
105 | char* marshalled_dataset=NULL;
|
---|
106 | int i;
|
---|
107 |
|
---|
108 | /*recover marshalled_dataset: */
|
---|
109 | marshalled_dataset=*pmarshalled_dataset;
|
---|
110 |
|
---|
111 | /*this time, no need to get enum type, the pointer directly points to the beginning of the
|
---|
112 | *object data (thanks to DataSet::Demarshall):*/
|
---|
113 | memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
|
---|
114 | memcpy(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);
|
---|
115 |
|
---|
116 | /*return: */
|
---|
117 | *pmarshalled_dataset=marshalled_dataset;
|
---|
118 | return;
|
---|
119 | }
|
---|
120 | /*}}}*/
|
---|
121 | /*FUNCTION PentaVertexInput::Enum{{{1*/
|
---|
122 | int PentaVertexInput::Enum(void){
|
---|
123 |
|
---|
124 | return PentaVertexInputEnum;
|
---|
125 |
|
---|
126 | }
|
---|
127 | /*}}}*/
|
---|
128 |
|
---|
129 | /*PentaVertexInput management*/
|
---|
130 | /*FUNCTION PentaVertexInput::copy{{{1*/
|
---|
131 | Object* PentaVertexInput::copy() {
|
---|
132 |
|
---|
133 | return new PentaVertexInput(this->enum_type,this->values);
|
---|
134 |
|
---|
135 | }
|
---|
136 | /*}}}*/
|
---|
137 | /*FUNCTION PentaVertexInput::EnumType{{{1*/
|
---|
138 | int PentaVertexInput::EnumType(void){
|
---|
139 |
|
---|
140 | return this->enum_type;
|
---|
141 |
|
---|
142 | }
|
---|
143 | /*}}}*/
|
---|
144 | /*FUNCTION PentaVertexInput::SpawnBeamInput{{{1*/
|
---|
145 | Input* PentaVertexInput::SpawnBeamInput(int* indices){
|
---|
146 |
|
---|
147 | /*output*/
|
---|
148 | BeamVertexInput* outinput=NULL;
|
---|
149 | double newvalues[2];
|
---|
150 |
|
---|
151 | /*Loop over the new indices*/
|
---|
152 | for(int i=0;i<2;i++){
|
---|
153 |
|
---|
154 | /*Check index value*/
|
---|
155 | ISSMASSERT(indices[i]>=0 && indices[i]<6);
|
---|
156 |
|
---|
157 | /*Assign value to new input*/
|
---|
158 | newvalues[i]=this->values[indices[i]];
|
---|
159 | }
|
---|
160 |
|
---|
161 | /*Create new Beam input*/
|
---|
162 | outinput=new BeamVertexInput(this->enum_type,&newvalues[0]);
|
---|
163 |
|
---|
164 | /*Assign output*/
|
---|
165 | return outinput;
|
---|
166 |
|
---|
167 | }
|
---|
168 | /*}}}*/
|
---|
169 | /*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
|
---|
170 | Input* PentaVertexInput::SpawnTriaInput(int* indices){
|
---|
171 |
|
---|
172 | /*output*/
|
---|
173 | TriaVertexInput* outinput=NULL;
|
---|
174 | double newvalues[3];
|
---|
175 |
|
---|
176 | /*Loop over the new indices*/
|
---|
177 | for(int i=0;i<3;i++){
|
---|
178 |
|
---|
179 | /*Check index value*/
|
---|
180 | ISSMASSERT(indices[i]>=0 && indices[i]<6);
|
---|
181 |
|
---|
182 | /*Assign value to new input*/
|
---|
183 | newvalues[i]=this->values[indices[i]];
|
---|
184 | }
|
---|
185 |
|
---|
186 | /*Create new Tria input*/
|
---|
187 | outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);
|
---|
188 |
|
---|
189 | /*Assign output*/
|
---|
190 | return outinput;
|
---|
191 |
|
---|
192 | }
|
---|
193 | /*}}}*/
|
---|
194 | /*FUNCTION PentaVertexInput::SpawnResult{{{1*/
|
---|
195 | ElementResult* PentaVertexInput::SpawnResult(int step, double time){
|
---|
196 |
|
---|
197 | return new PentaVertexElementResult(this->enum_type,this->values,step,time);
|
---|
198 |
|
---|
199 | }
|
---|
200 | /*}}}*/
|
---|
201 |
|
---|
202 | /*Object functions*/
|
---|
203 | /*FUNCTION PentaVertexInput::GetParameterValue(bool* pvalue) {{{1*/
|
---|
204 | void PentaVertexInput::GetParameterValue(bool* pvalue){ISSMERROR(" not supported yet!");}
|
---|
205 | /*}}}*/
|
---|
206 | /*FUNCTION PentaVertexInput::GetParameterValue(int* pvalue){{{1*/
|
---|
207 | void PentaVertexInput::GetParameterValue(int* pvalue){ISSMERROR(" not supported yet!");}
|
---|
208 | /*}}}*/
|
---|
209 | /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue){{{1*/
|
---|
210 | void PentaVertexInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
|
---|
211 | /*}}}*/
|
---|
212 | /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
|
---|
213 | void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
|
---|
214 | /*P1 interpolation on Gauss point*/
|
---|
215 |
|
---|
216 | /*intermediary*/
|
---|
217 | double l1l6[6];
|
---|
218 |
|
---|
219 | /*nodal functions: */
|
---|
220 | GetNodalFunctionsP1(&l1l6[0],gauss);
|
---|
221 |
|
---|
222 | /*Assign output pointers:*/
|
---|
223 | *pvalue=l1l6[0]*values[0]+l1l6[1]*values[1]+l1l6[2]*values[2]+l1l6[3]*values[3]+l1l6[4]*values[4]+l1l6[5]*values[5];
|
---|
224 |
|
---|
225 | }
|
---|
226 | /*}}}*/
|
---|
227 | /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
|
---|
228 | void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){ISSMERROR(" not supported yet!");}
|
---|
229 | /*}}}*/
|
---|
230 | /*FUNCTION PentaVertexInput::GetParameterValues{{{1*/
|
---|
231 | void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
|
---|
232 | /*It is assumed that output values has been correctly allocated*/
|
---|
233 |
|
---|
234 | int i,j;
|
---|
235 | double gauss[4];
|
---|
236 |
|
---|
237 | for (i=0;i<numgauss;i++){
|
---|
238 |
|
---|
239 | /*Get current Gauss point coordinates*/
|
---|
240 | for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];
|
---|
241 |
|
---|
242 | /*Assign parameter value*/
|
---|
243 | GetParameterValue(&values[i],&gauss[0]);
|
---|
244 | }
|
---|
245 | }
|
---|
246 | /*}}}*/
|
---|
247 | /*FUNCTION PentaVertexInput::GetParameterDerivativeValue{{{1*/
|
---|
248 | void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
|
---|
249 | /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
|
---|
250 | * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
|
---|
251 | * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
|
---|
252 | * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
|
---|
253 | *
|
---|
254 | * p is a vector of size 3x1 already allocated.
|
---|
255 | */
|
---|
256 |
|
---|
257 | const int NDOF3=3;
|
---|
258 | const int numgrids=6;
|
---|
259 | double dh1dh6[NDOF3][numgrids];
|
---|
260 |
|
---|
261 | /*Get nodal funnctions derivatives in actual coordinate system: */
|
---|
262 | GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
|
---|
263 |
|
---|
264 | p[0]=this->values[0]*dh1dh6[0][0]+this->values[1]*dh1dh6[0][1]+this->values[2]*dh1dh6[0][2]+this->values[3]*dh1dh6[0][3]+this->values[4]*dh1dh6[0][4]+this->values[5]*dh1dh6[0][5];
|
---|
265 | p[1]=this->values[0]*dh1dh6[1][0]+this->values[1]*dh1dh6[1][1]+this->values[2]*dh1dh6[1][2]+this->values[3]*dh1dh6[1][3]+this->values[4]*dh1dh6[1][4]+this->values[5]*dh1dh6[1][5];
|
---|
266 | p[2]=this->values[0]*dh1dh6[2][0]+this->values[1]*dh1dh6[2][1]+this->values[2]*dh1dh6[2][2]+this->values[3]*dh1dh6[2][3]+this->values[4]*dh1dh6[2][4]+this->values[5]*dh1dh6[2][5];
|
---|
267 |
|
---|
268 | }
|
---|
269 | /*}}}*/
|
---|
270 | /*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
|
---|
271 | void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
|
---|
272 | int i,j;
|
---|
273 |
|
---|
274 | const int numgrids=6;
|
---|
275 | const int DOFVELOCITY=3;
|
---|
276 | double B[8][27];
|
---|
277 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
278 | double velocity[numgrids][DOFVELOCITY];
|
---|
279 |
|
---|
280 | /*Get B matrix: */
|
---|
281 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
282 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
283 | for (i=0;i<6;i++){
|
---|
284 | for (j=0;j<3;j++){
|
---|
285 | B_reduced[i][j]=B[i][j];
|
---|
286 | }
|
---|
287 | for (j=4;j<7;j++){
|
---|
288 | B_reduced[i][j-1]=B[i][j];
|
---|
289 | }
|
---|
290 | for (j=8;j<11;j++){
|
---|
291 | B_reduced[i][j-2]=B[i][j];
|
---|
292 | }
|
---|
293 | for (j=12;j<15;j++){
|
---|
294 | B_reduced[i][j-3]=B[i][j];
|
---|
295 | }
|
---|
296 | for (j=16;j<19;j++){
|
---|
297 | B_reduced[i][j-4]=B[i][j];
|
---|
298 | }
|
---|
299 | for (j=20;j<23;j++){
|
---|
300 | B_reduced[i][j-5]=B[i][j];
|
---|
301 | }
|
---|
302 | }
|
---|
303 |
|
---|
304 | /*Here, we are computing the strain rate of (vx,0,0)*/
|
---|
305 | for(i=0;i<numgrids;i++){
|
---|
306 | velocity[i][0]=this->values[i];
|
---|
307 | velocity[i][1]=0.0;
|
---|
308 | velocity[i][2]=0.0;
|
---|
309 | }
|
---|
310 | /*Multiply B by velocity, to get strain rate: */
|
---|
311 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
|
---|
312 |
|
---|
313 | }
|
---|
314 | /*}}}*/
|
---|
315 | /*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
|
---|
316 | void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
|
---|
317 | int i,j;
|
---|
318 |
|
---|
319 | const int numgrids=6;
|
---|
320 | const int DOFVELOCITY=3;
|
---|
321 | double B[8][27];
|
---|
322 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
323 | double velocity[numgrids][DOFVELOCITY];
|
---|
324 |
|
---|
325 | /*Get B matrix: */
|
---|
326 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
327 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
328 | for (i=0;i<6;i++){
|
---|
329 | for (j=0;j<3;j++){
|
---|
330 | B_reduced[i][j]=B[i][j];
|
---|
331 | }
|
---|
332 | for (j=4;j<7;j++){
|
---|
333 | B_reduced[i][j-1]=B[i][j];
|
---|
334 | }
|
---|
335 | for (j=8;j<11;j++){
|
---|
336 | B_reduced[i][j-2]=B[i][j];
|
---|
337 | }
|
---|
338 | for (j=12;j<15;j++){
|
---|
339 | B_reduced[i][j-3]=B[i][j];
|
---|
340 | }
|
---|
341 | for (j=16;j<19;j++){
|
---|
342 | B_reduced[i][j-4]=B[i][j];
|
---|
343 | }
|
---|
344 | for (j=20;j<23;j++){
|
---|
345 | B_reduced[i][j-5]=B[i][j];
|
---|
346 | }
|
---|
347 | }
|
---|
348 |
|
---|
349 | /*Here, we are computing the strain rate of (0,vy,0)*/
|
---|
350 | for(i=0;i<numgrids;i++){
|
---|
351 | velocity[i][0]=0.0;
|
---|
352 | velocity[i][1]=this->values[i];
|
---|
353 | velocity[i][2]=0.0;
|
---|
354 | }
|
---|
355 | /*Multiply B by velocity, to get strain rate: */
|
---|
356 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
|
---|
357 |
|
---|
358 | }
|
---|
359 | /*}}}*/
|
---|
360 | /*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
|
---|
361 | void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
|
---|
362 | int i,j;
|
---|
363 |
|
---|
364 | const int numgrids=6;
|
---|
365 | const int DOFVELOCITY=3;
|
---|
366 | double B[8][27];
|
---|
367 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
368 | double velocity[numgrids][DOFVELOCITY];
|
---|
369 |
|
---|
370 | /*Get B matrix: */
|
---|
371 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
372 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
373 | for (i=0;i<6;i++){
|
---|
374 | for (j=0;j<3;j++){
|
---|
375 | B_reduced[i][j]=B[i][j];
|
---|
376 | }
|
---|
377 | for (j=4;j<7;j++){
|
---|
378 | B_reduced[i][j-1]=B[i][j];
|
---|
379 | }
|
---|
380 | for (j=8;j<11;j++){
|
---|
381 | B_reduced[i][j-2]=B[i][j];
|
---|
382 | }
|
---|
383 | for (j=12;j<15;j++){
|
---|
384 | B_reduced[i][j-3]=B[i][j];
|
---|
385 | }
|
---|
386 | for (j=16;j<19;j++){
|
---|
387 | B_reduced[i][j-4]=B[i][j];
|
---|
388 | }
|
---|
389 | for (j=20;j<23;j++){
|
---|
390 | B_reduced[i][j-5]=B[i][j];
|
---|
391 | }
|
---|
392 | }
|
---|
393 |
|
---|
394 | /*Here, we are computing the strain rate of (0,0,vz)*/
|
---|
395 | for(i=0;i<numgrids;i++){
|
---|
396 | velocity[i][0]=0.0;
|
---|
397 | velocity[i][1]=0.0;
|
---|
398 | velocity[i][2]=this->values[i];
|
---|
399 | }
|
---|
400 |
|
---|
401 | /*Multiply B by velocity, to get strain rate: */
|
---|
402 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
|
---|
403 |
|
---|
404 | }
|
---|
405 | /*}}}*/
|
---|
406 | /*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
|
---|
407 | void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
|
---|
408 |
|
---|
409 | int i;
|
---|
410 | const int numgrids=6;
|
---|
411 | const int NDOF2=2;
|
---|
412 | double B[5][NDOF2*numgrids];
|
---|
413 | double velocity[numgrids][NDOF2];
|
---|
414 |
|
---|
415 | /*Get B matrix: */
|
---|
416 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
417 |
|
---|
418 | /*Here, we are computing the strain rate of (vx,0)*/
|
---|
419 | for(i=0;i<numgrids;i++){
|
---|
420 | velocity[i][0]=this->values[i];
|
---|
421 | velocity[i][1]=0.0;
|
---|
422 | }
|
---|
423 |
|
---|
424 | /*Multiply B by velocity, to get strain rate: */
|
---|
425 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
426 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
427 | epsilonvx,0);
|
---|
428 |
|
---|
429 | }
|
---|
430 | /*}}}*/
|
---|
431 | /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
|
---|
432 | void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
|
---|
433 |
|
---|
434 | int i;
|
---|
435 | const int numgrids=6;
|
---|
436 | const int NDOF2=2;
|
---|
437 | double B[5][NDOF2*numgrids];
|
---|
438 | double velocity[numgrids][NDOF2];
|
---|
439 |
|
---|
440 | /*Get B matrix: */
|
---|
441 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
442 |
|
---|
443 | /*Here, we are computing the strain rate of (0,vy)*/
|
---|
444 | for(i=0;i<numgrids;i++){
|
---|
445 | velocity[i][0]=0.0;
|
---|
446 | velocity[i][1]=this->values[i];
|
---|
447 | }
|
---|
448 |
|
---|
449 | /*Multiply B by velocity, to get strain rate: */
|
---|
450 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
451 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
452 | epsilonvy,0);
|
---|
453 |
|
---|
454 | }
|
---|
455 | /*}}}*/
|
---|
456 | /*FUNCTION PentaVertexInput::ChangeEnum{{{1*/
|
---|
457 | void PentaVertexInput::ChangeEnum(int newenumtype){
|
---|
458 | this->enum_type=newenumtype;
|
---|
459 | }
|
---|
460 | /*}}}*/
|
---|
461 | /*FUNCTION PentaVertexInput::GetParameterAverage{{{1*/
|
---|
462 | void PentaVertexInput::GetParameterAverage(double* pvalue){
|
---|
463 | *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
|
---|
464 | }
|
---|
465 | /*}}}*/
|
---|
466 |
|
---|
467 | /*Intermediary*/
|
---|
468 | /*FUNCTION PentaVertexInput::SquareMin{{{1*/
|
---|
469 | void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
|
---|
470 |
|
---|
471 | int i;
|
---|
472 | const int numnodes=6;
|
---|
473 | double valuescopy[numnodes];
|
---|
474 | double squaremin;
|
---|
475 |
|
---|
476 | /*First, copy values, to process units if requested: */
|
---|
477 | for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
|
---|
478 |
|
---|
479 | /*Process units if requested: */
|
---|
480 | if(process_units)NodalValuesUnitConversion(&valuescopy[0],numnodes,enum_type,parameters);
|
---|
481 |
|
---|
482 | /*Now, figure out minimum of valuescopy: */
|
---|
483 | squaremin=pow(valuescopy[0],2);
|
---|
484 | for(i=1;i<numnodes;i++){
|
---|
485 | if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
|
---|
486 | }
|
---|
487 | /*Assign output pointers:*/
|
---|
488 | *psquaremin=squaremin;
|
---|
489 | }
|
---|
490 | /*}}}*/
|
---|
491 | /*FUNCTION PentaVertexInput::Scale{{{1*/
|
---|
492 | void PentaVertexInput::Scale(double scale_factor){
|
---|
493 |
|
---|
494 | int i;
|
---|
495 | const int numgrids=6;
|
---|
496 |
|
---|
497 | for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
|
---|
498 | }
|
---|
499 | /*}}}*/
|
---|
500 | /*FUNCTION PentaVertexInput::AXPY{{{1*/
|
---|
501 | void PentaVertexInput::AXPY(Input* xinput,double scalar){
|
---|
502 |
|
---|
503 | int i;
|
---|
504 | const int numgrids=6;
|
---|
505 | PentaVertexInput* xpentavertexinput=NULL;
|
---|
506 |
|
---|
507 | /*xinput is of the same type, so cast it: */
|
---|
508 | xpentavertexinput=(PentaVertexInput*)xinput;
|
---|
509 |
|
---|
510 | /*Carry out the AXPY operation depending on type:*/
|
---|
511 | switch(xinput->Enum()){
|
---|
512 |
|
---|
513 | case PentaVertexInputEnum:
|
---|
514 | for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
|
---|
515 | return;
|
---|
516 |
|
---|
517 | default:
|
---|
518 | ISSMERROR("not implemented yet");
|
---|
519 | }
|
---|
520 |
|
---|
521 | }
|
---|
522 | /*}}}*/
|
---|
523 | /*FUNCTION PentaVertexInput::Constrain{{{1*/
|
---|
524 | void PentaVertexInput::Constrain(double cm_min, double cm_max){
|
---|
525 |
|
---|
526 | int i;
|
---|
527 | const int numgrids=6;
|
---|
528 |
|
---|
529 | if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
|
---|
530 | if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
|
---|
531 |
|
---|
532 | }
|
---|
533 | /*}}}*/
|
---|
534 | /*FUNCTION PentaVertexInput::Extrude{{{1*/
|
---|
535 | void PentaVertexInput::Extrude(void){
|
---|
536 |
|
---|
537 | int i;
|
---|
538 |
|
---|
539 | /*First 3 values copied on 3 last values*/
|
---|
540 | for(i=0;i<3;i++) this->values[3+i]=this->values[i];
|
---|
541 |
|
---|
542 | }
|
---|
543 | /*}}}*/
|
---|
544 | /*FUNCTION PentaVertexInput::VerticallyIntegrate{{{1*/
|
---|
545 | void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){
|
---|
546 |
|
---|
547 | /*Intermediaries*/
|
---|
548 | int i;
|
---|
549 | const int numgrids = 6;
|
---|
550 | int num_thickness_values;
|
---|
551 | double *thickness_values = NULL;
|
---|
552 |
|
---|
553 | /*Check that input provided is a thickness*/
|
---|
554 | if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumAsString(thickness_input->EnumType()));
|
---|
555 |
|
---|
556 | /*Get Thickness value pointer*/
|
---|
557 | thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
|
---|
558 |
|
---|
559 | /*vertically integrate depending on type:*/
|
---|
560 | switch(thickness_input->Enum()){
|
---|
561 |
|
---|
562 | case PentaVertexInputEnum:
|
---|
563 | for(i=0;i<3;i++){
|
---|
564 | this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
|
---|
565 | this->values[i+3]=this->values[i];
|
---|
566 | }
|
---|
567 | return;
|
---|
568 |
|
---|
569 | default:
|
---|
570 | ISSMERROR("not implemented yet");
|
---|
571 | }
|
---|
572 | }
|
---|
573 | /*}}}*/
|
---|
574 | /*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
|
---|
575 | Input* PentaVertexInput::PointwiseDivide(Input* inputB){
|
---|
576 |
|
---|
577 | /*Ouput*/
|
---|
578 | PentaVertexInput* outinput=NULL;
|
---|
579 |
|
---|
580 | /*Intermediaries*/
|
---|
581 | int i;
|
---|
582 | PentaVertexInput *xinputB = NULL;
|
---|
583 | int B_numvalues;
|
---|
584 | double *B_values = NULL;
|
---|
585 | const int numgrids = 6;
|
---|
586 | double AdotBvalues[numgrids];
|
---|
587 |
|
---|
588 | /*Check that inputB is of the same type*/
|
---|
589 | if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumAsString(inputB->Enum()));
|
---|
590 | xinputB=(PentaVertexInput*)inputB;
|
---|
591 |
|
---|
592 | /*Create point wise sum*/
|
---|
593 | for(i=0;i<numgrids;i++){
|
---|
594 | ISSMASSERT(xinputB->values[i]!=0);
|
---|
595 | AdotBvalues[i]=this->values[i]/xinputB->values[i];
|
---|
596 | }
|
---|
597 |
|
---|
598 | /*Create new Penta vertex input (copy of current input)*/
|
---|
599 | outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]);
|
---|
600 |
|
---|
601 | /*Return output pointer*/
|
---|
602 | return outinput;
|
---|
603 |
|
---|
604 | }
|
---|
605 | /*}}}*/
|
---|
606 | /*FUNCTION PentaVertexInput::GetVectorFromInputs{{{1*/
|
---|
607 | void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
|
---|
608 |
|
---|
609 | const int numvertices=6;
|
---|
610 | VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);
|
---|
611 |
|
---|
612 | } /*}}}*/
|
---|
613 | /*FUNCTION PentaVertexInput::GetValuesPtr{{{1*/
|
---|
614 | void PentaVertexInput::GetValuesPtr(double** pvalues,int* pnum_values){
|
---|
615 |
|
---|
616 | *pvalues=this->values;
|
---|
617 | *pnum_values=6;
|
---|
618 |
|
---|
619 | }
|
---|
620 | /*}}}*/
|
---|