[3683] | 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"
|
---|
[4236] | 16 | #include "../../Container/Container.h"
|
---|
[3775] | 17 | #include "../../include/include.h"
|
---|
[3683] | 18 |
|
---|
[4248] | 19 | /*PentaVertexInput constructors and destructor*/
|
---|
[3683] | 20 | /*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
|
---|
| 21 | PentaVertexInput::PentaVertexInput(){
|
---|
| 22 | return;
|
---|
| 23 | }
|
---|
| 24 | /*}}}*/
|
---|
[3847] | 25 | /*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
|
---|
[4882] | 26 | PentaVertexInput::PentaVertexInput(int in_enum_type,double* in_values)
|
---|
| 27 | :PentaRef(1)
|
---|
| 28 | {
|
---|
[3683] | 29 |
|
---|
[4882] | 30 | /*Set PentaRef*/
|
---|
| 31 | this->SetElementType(P1Enum,0);
|
---|
| 32 | this->element_type=P1Enum;
|
---|
| 33 |
|
---|
[3683] | 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 |
|
---|
[4248] | 49 | /*Object virtual functions definitions:*/
|
---|
| 50 | /*FUNCTION PentaVertexInput::Echo {{{1*/
|
---|
| 51 | void PentaVertexInput::Echo(void){
|
---|
| 52 | this->DeepEcho();
|
---|
[3683] | 53 | }
|
---|
| 54 | /*}}}*/
|
---|
| 55 | /*FUNCTION PentaVertexInput::DeepEcho{{{1*/
|
---|
| 56 | void PentaVertexInput::DeepEcho(void){
|
---|
| 57 |
|
---|
| 58 | printf("PentaVertexInput:\n");
|
---|
[5103] | 59 | printf(" enum: %i (%s)\n",this->enum_type,EnumToString(this->enum_type));
|
---|
[3847] | 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]);
|
---|
[3683] | 61 | }
|
---|
| 62 | /*}}}*/
|
---|
[4248] | 63 | /*FUNCTION PentaVertexInput::Id{{{1*/
|
---|
| 64 | int PentaVertexInput::Id(void){ return -1; }
|
---|
[3683] | 65 | /*}}}*/
|
---|
[4248] | 66 | /*FUNCTION PentaVertexInput::MyRank{{{1*/
|
---|
| 67 | int PentaVertexInput::MyRank(void){
|
---|
| 68 | extern int my_rank;
|
---|
| 69 | return my_rank;
|
---|
[3683] | 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 | /*}}}*/
|
---|
[4248] | 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;
|
---|
[3683] | 119 | }
|
---|
| 120 | /*}}}*/
|
---|
[4248] | 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 | /*}}}*/
|
---|
[3847] | 144 | /*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
|
---|
| 145 | Input* PentaVertexInput::SpawnTriaInput(int* indices){
|
---|
[3683] | 146 |
|
---|
[3847] | 147 | /*output*/
|
---|
| 148 | TriaVertexInput* outinput=NULL;
|
---|
| 149 | double newvalues[3];
|
---|
| 150 |
|
---|
| 151 | /*Loop over the new indices*/
|
---|
| 152 | for(int i=0;i<3;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 Tria input*/
|
---|
| 162 | outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);
|
---|
| 163 |
|
---|
| 164 | /*Assign output*/
|
---|
| 165 | return outinput;
|
---|
| 166 |
|
---|
| 167 | }
|
---|
| 168 | /*}}}*/
|
---|
[4037] | 169 | /*FUNCTION PentaVertexInput::SpawnResult{{{1*/
|
---|
[4050] | 170 | ElementResult* PentaVertexInput::SpawnResult(int step, double time){
|
---|
[3847] | 171 |
|
---|
[4050] | 172 | return new PentaVertexElementResult(this->enum_type,this->values,step,time);
|
---|
[4037] | 173 |
|
---|
| 174 | }
|
---|
| 175 | /*}}}*/
|
---|
| 176 |
|
---|
[3683] | 177 | /*Object functions*/
|
---|
[5629] | 178 | /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
|
---|
[3840] | 179 | void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
|
---|
| 180 |
|
---|
[4921] | 181 | /*Call PentaRef function*/
|
---|
| 182 | PentaRef::GetParameterValue(pvalue,&values[0],gauss);
|
---|
[3840] | 183 |
|
---|
| 184 | }
|
---|
[3683] | 185 | /*}}}*/
|
---|
[5647] | 186 | /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
|
---|
| 187 | void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){
|
---|
[5629] | 188 |
|
---|
| 189 | /*Call PentaRef function*/
|
---|
| 190 | PentaRef::GetParameterValue(pvalue,&values[0],gauss);
|
---|
| 191 |
|
---|
| 192 | }
|
---|
| 193 | /*}}}*/
|
---|
[4546] | 194 | /*FUNCTION PentaVertexInput::GetParameterValues{{{1*/
|
---|
[3840] | 195 | void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
|
---|
| 196 | /*It is assumed that output values has been correctly allocated*/
|
---|
| 197 |
|
---|
| 198 | int i,j;
|
---|
| 199 | double gauss[4];
|
---|
| 200 |
|
---|
| 201 | for (i=0;i<numgauss;i++){
|
---|
| 202 |
|
---|
| 203 | /*Get current Gauss point coordinates*/
|
---|
| 204 | for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];
|
---|
| 205 |
|
---|
| 206 | /*Assign parameter value*/
|
---|
| 207 | GetParameterValue(&values[i],&gauss[0]);
|
---|
| 208 | }
|
---|
| 209 | }
|
---|
[3683] | 210 | /*}}}*/
|
---|
[5629] | 211 | /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
|
---|
[3840] | 212 | void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
|
---|
| 213 |
|
---|
[4921] | 214 | /*Call PentaRef function*/
|
---|
| 215 | PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
|
---|
[3840] | 216 | }
|
---|
[3683] | 217 | /*}}}*/
|
---|
[5647] | 218 | /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/
|
---|
| 219 | void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){
|
---|
[5629] | 220 |
|
---|
| 221 | /*Call PentaRef function*/
|
---|
| 222 | PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
|
---|
| 223 | }
|
---|
| 224 | /*}}}*/
|
---|
[4546] | 225 | /*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
|
---|
[3855] | 226 | void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
|
---|
[3840] | 227 | int i,j;
|
---|
| 228 |
|
---|
| 229 | const int numgrids=6;
|
---|
| 230 | const int DOFVELOCITY=3;
|
---|
| 231 | double B[8][27];
|
---|
| 232 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
[3875] | 233 | double velocity[numgrids][DOFVELOCITY];
|
---|
[3840] | 234 |
|
---|
| 235 | /*Get B matrix: */
|
---|
| 236 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
| 237 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
| 238 | for (i=0;i<6;i++){
|
---|
| 239 | for (j=0;j<3;j++){
|
---|
| 240 | B_reduced[i][j]=B[i][j];
|
---|
| 241 | }
|
---|
| 242 | for (j=4;j<7;j++){
|
---|
| 243 | B_reduced[i][j-1]=B[i][j];
|
---|
| 244 | }
|
---|
| 245 | for (j=8;j<11;j++){
|
---|
| 246 | B_reduced[i][j-2]=B[i][j];
|
---|
| 247 | }
|
---|
| 248 | for (j=12;j<15;j++){
|
---|
| 249 | B_reduced[i][j-3]=B[i][j];
|
---|
| 250 | }
|
---|
| 251 | for (j=16;j<19;j++){
|
---|
| 252 | B_reduced[i][j-4]=B[i][j];
|
---|
| 253 | }
|
---|
| 254 | for (j=20;j<23;j++){
|
---|
| 255 | B_reduced[i][j-5]=B[i][j];
|
---|
| 256 | }
|
---|
| 257 | }
|
---|
| 258 |
|
---|
| 259 | /*Here, we are computing the strain rate of (vx,0,0)*/
|
---|
| 260 | for(i=0;i<numgrids;i++){
|
---|
| 261 | velocity[i][0]=this->values[i];
|
---|
| 262 | velocity[i][1]=0.0;
|
---|
| 263 | velocity[i][2]=0.0;
|
---|
| 264 | }
|
---|
| 265 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 266 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
|
---|
| 267 |
|
---|
| 268 | }
|
---|
| 269 | /*}}}*/
|
---|
[4546] | 270 | /*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
|
---|
[3855] | 271 | void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
|
---|
[3840] | 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];
|
---|
[3875] | 278 | double velocity[numgrids][DOFVELOCITY];
|
---|
[3840] | 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 (0,vy,0)*/
|
---|
| 305 | for(i=0;i<numgrids;i++){
|
---|
| 306 | velocity[i][0]=0.0;
|
---|
| 307 | velocity[i][1]=this->values[i];
|
---|
| 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,epsilonvy,0);
|
---|
| 312 |
|
---|
| 313 | }
|
---|
| 314 | /*}}}*/
|
---|
[4546] | 315 | /*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
|
---|
[3855] | 316 | void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
|
---|
[3840] | 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];
|
---|
[3875] | 323 | double velocity[numgrids][DOFVELOCITY];
|
---|
[3840] | 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,0,vz)*/
|
---|
| 350 | for(i=0;i<numgrids;i++){
|
---|
| 351 | velocity[i][0]=0.0;
|
---|
| 352 | velocity[i][1]=0.0;
|
---|
| 353 | velocity[i][2]=this->values[i];
|
---|
| 354 | }
|
---|
| 355 |
|
---|
| 356 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 357 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
|
---|
| 358 |
|
---|
| 359 | }
|
---|
| 360 | /*}}}*/
|
---|
[4546] | 361 | /*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
|
---|
[3855] | 362 | void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){
|
---|
[3840] | 363 |
|
---|
[3855] | 364 | int i;
|
---|
| 365 | const int numgrids=6;
|
---|
| 366 | const int NDOF2=2;
|
---|
| 367 | double B[5][NDOF2*numgrids];
|
---|
| 368 | double velocity[numgrids][NDOF2];
|
---|
| 369 |
|
---|
| 370 | /*Get B matrix: */
|
---|
| 371 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
| 372 |
|
---|
| 373 | /*Here, we are computing the strain rate of (vx,0)*/
|
---|
| 374 | for(i=0;i<numgrids;i++){
|
---|
| 375 | velocity[i][0]=this->values[i];
|
---|
| 376 | velocity[i][1]=0.0;
|
---|
[3840] | 377 | }
|
---|
| 378 |
|
---|
[3855] | 379 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 380 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
| 381 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
| 382 | epsilonvx,0);
|
---|
| 383 |
|
---|
[3840] | 384 | }
|
---|
| 385 | /*}}}*/
|
---|
[4546] | 386 | /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
|
---|
[3855] | 387 | void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
|
---|
| 388 |
|
---|
| 389 | int i;
|
---|
| 390 | const int numgrids=6;
|
---|
| 391 | const int NDOF2=2;
|
---|
| 392 | double B[5][NDOF2*numgrids];
|
---|
| 393 | double velocity[numgrids][NDOF2];
|
---|
| 394 |
|
---|
| 395 | /*Get B matrix: */
|
---|
| 396 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
| 397 |
|
---|
| 398 | /*Here, we are computing the strain rate of (0,vy)*/
|
---|
| 399 | for(i=0;i<numgrids;i++){
|
---|
| 400 | velocity[i][0]=0.0;
|
---|
| 401 | velocity[i][1]=this->values[i];
|
---|
| 402 | }
|
---|
| 403 |
|
---|
| 404 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 405 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
| 406 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
| 407 | epsilonvy,0);
|
---|
| 408 |
|
---|
| 409 | }
|
---|
| 410 | /*}}}*/
|
---|
[5647] | 411 | /*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
|
---|
| 412 | void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
|
---|
| 413 | int i,j;
|
---|
| 414 |
|
---|
| 415 | const int numgrids=6;
|
---|
| 416 | const int DOFVELOCITY=3;
|
---|
| 417 | double B[8][27];
|
---|
| 418 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
| 419 | double velocity[numgrids][DOFVELOCITY];
|
---|
| 420 |
|
---|
| 421 | /*Get B matrix: */
|
---|
| 422 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
| 423 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
| 424 | for (i=0;i<6;i++){
|
---|
| 425 | for (j=0;j<3;j++){
|
---|
| 426 | B_reduced[i][j]=B[i][j];
|
---|
| 427 | }
|
---|
| 428 | for (j=4;j<7;j++){
|
---|
| 429 | B_reduced[i][j-1]=B[i][j];
|
---|
| 430 | }
|
---|
| 431 | for (j=8;j<11;j++){
|
---|
| 432 | B_reduced[i][j-2]=B[i][j];
|
---|
| 433 | }
|
---|
| 434 | for (j=12;j<15;j++){
|
---|
| 435 | B_reduced[i][j-3]=B[i][j];
|
---|
| 436 | }
|
---|
| 437 | for (j=16;j<19;j++){
|
---|
| 438 | B_reduced[i][j-4]=B[i][j];
|
---|
| 439 | }
|
---|
| 440 | for (j=20;j<23;j++){
|
---|
| 441 | B_reduced[i][j-5]=B[i][j];
|
---|
| 442 | }
|
---|
| 443 | }
|
---|
| 444 |
|
---|
| 445 | /*Here, we are computing the strain rate of (vx,0,0)*/
|
---|
| 446 | for(i=0;i<numgrids;i++){
|
---|
| 447 | velocity[i][0]=this->values[i];
|
---|
| 448 | velocity[i][1]=0.0;
|
---|
| 449 | velocity[i][2]=0.0;
|
---|
| 450 | }
|
---|
| 451 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 452 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
|
---|
| 453 |
|
---|
| 454 | }
|
---|
| 455 | /*}}}*/
|
---|
| 456 | /*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
|
---|
| 457 | void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
|
---|
| 458 | int i,j;
|
---|
| 459 |
|
---|
| 460 | const int numgrids=6;
|
---|
| 461 | const int DOFVELOCITY=3;
|
---|
| 462 | double B[8][27];
|
---|
| 463 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
| 464 | double velocity[numgrids][DOFVELOCITY];
|
---|
| 465 |
|
---|
| 466 | /*Get B matrix: */
|
---|
| 467 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
| 468 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
| 469 | for (i=0;i<6;i++){
|
---|
| 470 | for (j=0;j<3;j++){
|
---|
| 471 | B_reduced[i][j]=B[i][j];
|
---|
| 472 | }
|
---|
| 473 | for (j=4;j<7;j++){
|
---|
| 474 | B_reduced[i][j-1]=B[i][j];
|
---|
| 475 | }
|
---|
| 476 | for (j=8;j<11;j++){
|
---|
| 477 | B_reduced[i][j-2]=B[i][j];
|
---|
| 478 | }
|
---|
| 479 | for (j=12;j<15;j++){
|
---|
| 480 | B_reduced[i][j-3]=B[i][j];
|
---|
| 481 | }
|
---|
| 482 | for (j=16;j<19;j++){
|
---|
| 483 | B_reduced[i][j-4]=B[i][j];
|
---|
| 484 | }
|
---|
| 485 | for (j=20;j<23;j++){
|
---|
| 486 | B_reduced[i][j-5]=B[i][j];
|
---|
| 487 | }
|
---|
| 488 | }
|
---|
| 489 |
|
---|
| 490 | /*Here, we are computing the strain rate of (0,vy,0)*/
|
---|
| 491 | for(i=0;i<numgrids;i++){
|
---|
| 492 | velocity[i][0]=0.0;
|
---|
| 493 | velocity[i][1]=this->values[i];
|
---|
| 494 | velocity[i][2]=0.0;
|
---|
| 495 | }
|
---|
| 496 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 497 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
|
---|
| 498 |
|
---|
| 499 | }
|
---|
| 500 | /*}}}*/
|
---|
| 501 | /*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
|
---|
| 502 | void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){
|
---|
| 503 | int i,j;
|
---|
| 504 |
|
---|
| 505 | const int numgrids=6;
|
---|
| 506 | const int DOFVELOCITY=3;
|
---|
| 507 | double B[8][27];
|
---|
| 508 | double B_reduced[6][DOFVELOCITY*numgrids];
|
---|
| 509 | double velocity[numgrids][DOFVELOCITY];
|
---|
| 510 |
|
---|
| 511 | /*Get B matrix: */
|
---|
| 512 | GetBStokes(&B[0][0], xyz_list, gauss);
|
---|
| 513 | /*Create a reduced matrix of B to get rid of pressure */
|
---|
| 514 | for (i=0;i<6;i++){
|
---|
| 515 | for (j=0;j<3;j++){
|
---|
| 516 | B_reduced[i][j]=B[i][j];
|
---|
| 517 | }
|
---|
| 518 | for (j=4;j<7;j++){
|
---|
| 519 | B_reduced[i][j-1]=B[i][j];
|
---|
| 520 | }
|
---|
| 521 | for (j=8;j<11;j++){
|
---|
| 522 | B_reduced[i][j-2]=B[i][j];
|
---|
| 523 | }
|
---|
| 524 | for (j=12;j<15;j++){
|
---|
| 525 | B_reduced[i][j-3]=B[i][j];
|
---|
| 526 | }
|
---|
| 527 | for (j=16;j<19;j++){
|
---|
| 528 | B_reduced[i][j-4]=B[i][j];
|
---|
| 529 | }
|
---|
| 530 | for (j=20;j<23;j++){
|
---|
| 531 | B_reduced[i][j-5]=B[i][j];
|
---|
| 532 | }
|
---|
| 533 | }
|
---|
| 534 |
|
---|
| 535 | /*Here, we are computing the strain rate of (0,0,vz)*/
|
---|
| 536 | for(i=0;i<numgrids;i++){
|
---|
| 537 | velocity[i][0]=0.0;
|
---|
| 538 | velocity[i][1]=0.0;
|
---|
| 539 | velocity[i][2]=this->values[i];
|
---|
| 540 | }
|
---|
| 541 |
|
---|
| 542 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 543 | MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
|
---|
| 544 |
|
---|
| 545 | }
|
---|
| 546 | /*}}}*/
|
---|
| 547 | /*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
|
---|
| 548 | void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
|
---|
| 549 |
|
---|
| 550 | int i;
|
---|
| 551 | const int numgrids=6;
|
---|
| 552 | const int NDOF2=2;
|
---|
| 553 | double B[5][NDOF2*numgrids];
|
---|
| 554 | double velocity[numgrids][NDOF2];
|
---|
| 555 |
|
---|
| 556 | /*Get B matrix: */
|
---|
| 557 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
| 558 |
|
---|
| 559 | /*Here, we are computing the strain rate of (vx,0)*/
|
---|
| 560 | for(i=0;i<numgrids;i++){
|
---|
| 561 | velocity[i][0]=this->values[i];
|
---|
| 562 | velocity[i][1]=0.0;
|
---|
| 563 | }
|
---|
| 564 |
|
---|
| 565 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 566 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
| 567 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
| 568 | epsilonvx,0);
|
---|
| 569 |
|
---|
| 570 | }
|
---|
| 571 | /*}}}*/
|
---|
| 572 | /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
|
---|
| 573 | void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
|
---|
| 574 |
|
---|
| 575 | int i;
|
---|
| 576 | const int numgrids=6;
|
---|
| 577 | const int NDOF2=2;
|
---|
| 578 | double B[5][NDOF2*numgrids];
|
---|
| 579 | double velocity[numgrids][NDOF2];
|
---|
| 580 |
|
---|
| 581 | /*Get B matrix: */
|
---|
| 582 | GetBPattyn(&B[0][0], xyz_list, gauss);
|
---|
| 583 |
|
---|
| 584 | /*Here, we are computing the strain rate of (0,vy)*/
|
---|
| 585 | for(i=0;i<numgrids;i++){
|
---|
| 586 | velocity[i][0]=0.0;
|
---|
| 587 | velocity[i][1]=this->values[i];
|
---|
| 588 | }
|
---|
| 589 |
|
---|
| 590 | /*Multiply B by velocity, to get strain rate: */
|
---|
| 591 | MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
|
---|
| 592 | &velocity[0][0],NDOF2*numgrids,1,0,
|
---|
| 593 | epsilonvy,0);
|
---|
| 594 |
|
---|
| 595 | }
|
---|
| 596 | /*}}}*/
|
---|
[4546] | 597 | /*FUNCTION PentaVertexInput::ChangeEnum{{{1*/
|
---|
[3732] | 598 | void PentaVertexInput::ChangeEnum(int newenumtype){
|
---|
| 599 | this->enum_type=newenumtype;
|
---|
| 600 | }
|
---|
| 601 | /*}}}*/
|
---|
[4546] | 602 | /*FUNCTION PentaVertexInput::GetParameterAverage{{{1*/
|
---|
[3830] | 603 | void PentaVertexInput::GetParameterAverage(double* pvalue){
|
---|
| 604 | *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
|
---|
| 605 | }
|
---|
| 606 | /*}}}*/
|
---|
[3840] | 607 |
|
---|
| 608 | /*Intermediary*/
|
---|
[4471] | 609 | /*FUNCTION PentaVertexInput::SquareMin{{{1*/
|
---|
[4042] | 610 | void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
|
---|
| 611 |
|
---|
| 612 | int i;
|
---|
| 613 | const int numnodes=6;
|
---|
| 614 | double valuescopy[numnodes];
|
---|
| 615 | double squaremin;
|
---|
| 616 |
|
---|
| 617 | /*First, copy values, to process units if requested: */
|
---|
| 618 | for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
|
---|
| 619 |
|
---|
| 620 | /*Process units if requested: */
|
---|
[5529] | 621 | if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type,parameters);
|
---|
[4042] | 622 |
|
---|
| 623 | /*Now, figure out minimum of valuescopy: */
|
---|
| 624 | squaremin=pow(valuescopy[0],2);
|
---|
| 625 | for(i=1;i<numnodes;i++){
|
---|
| 626 | if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
|
---|
| 627 | }
|
---|
| 628 | /*Assign output pointers:*/
|
---|
| 629 | *psquaremin=squaremin;
|
---|
| 630 | }
|
---|
| 631 | /*}}}*/
|
---|
[5017] | 632 | /*FUNCTION PentaVertexInput::ConstrainMin{{{1*/
|
---|
| 633 | void PentaVertexInput::ConstrainMin(double minimum){
|
---|
| 634 |
|
---|
| 635 | int i;
|
---|
| 636 | const int numgrids=6;
|
---|
| 637 |
|
---|
| 638 | for(i=0;i<numgrids;i++) if (values[i]<minimum) values[i]=minimum;
|
---|
| 639 | }
|
---|
| 640 | /*}}}*/
|
---|
[5513] | 641 | /*FUNCTION PentaVertexInput::InfinityNorm{{{1*/
|
---|
| 642 | double PentaVertexInput::InfinityNorm(void){
|
---|
| 643 |
|
---|
| 644 | /*Output*/
|
---|
| 645 | const int numgrids=6;
|
---|
| 646 | double norm=0;
|
---|
| 647 |
|
---|
| 648 | for(int i=0;i<numgrids;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
|
---|
| 649 | return norm;
|
---|
| 650 | }
|
---|
| 651 | /*}}}*/
|
---|
[4471] | 652 | /*FUNCTION PentaVertexInput::Scale{{{1*/
|
---|
[4047] | 653 | void PentaVertexInput::Scale(double scale_factor){
|
---|
| 654 |
|
---|
| 655 | int i;
|
---|
| 656 | const int numgrids=6;
|
---|
| 657 |
|
---|
| 658 | for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
|
---|
| 659 | }
|
---|
| 660 | /*}}}*/
|
---|
[4471] | 661 | /*FUNCTION PentaVertexInput::AXPY{{{1*/
|
---|
[4048] | 662 | void PentaVertexInput::AXPY(Input* xinput,double scalar){
|
---|
| 663 |
|
---|
| 664 | int i;
|
---|
| 665 | const int numgrids=6;
|
---|
| 666 | PentaVertexInput* xpentavertexinput=NULL;
|
---|
| 667 |
|
---|
| 668 | /*xinput is of the same type, so cast it: */
|
---|
[4050] | 669 | xpentavertexinput=(PentaVertexInput*)xinput;
|
---|
[4048] | 670 |
|
---|
[4174] | 671 | /*Carry out the AXPY operation depending on type:*/
|
---|
| 672 | switch(xinput->Enum()){
|
---|
[4048] | 673 |
|
---|
[4174] | 674 | case PentaVertexInputEnum:
|
---|
| 675 | for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
|
---|
| 676 | return;
|
---|
| 677 |
|
---|
| 678 | default:
|
---|
| 679 | ISSMERROR("not implemented yet");
|
---|
| 680 | }
|
---|
| 681 |
|
---|
[4048] | 682 | }
|
---|
| 683 | /*}}}*/
|
---|
[4471] | 684 | /*FUNCTION PentaVertexInput::Constrain{{{1*/
|
---|
[4048] | 685 | void PentaVertexInput::Constrain(double cm_min, double cm_max){
|
---|
| 686 |
|
---|
| 687 | int i;
|
---|
| 688 | const int numgrids=6;
|
---|
| 689 |
|
---|
| 690 | if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
|
---|
| 691 | if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
|
---|
| 692 |
|
---|
| 693 | }
|
---|
| 694 | /*}}}*/
|
---|
[4471] | 695 | /*FUNCTION PentaVertexInput::Extrude{{{1*/
|
---|
[4274] | 696 | void PentaVertexInput::Extrude(void){
|
---|
| 697 |
|
---|
| 698 | int i;
|
---|
| 699 |
|
---|
| 700 | /*First 3 values copied on 3 last values*/
|
---|
| 701 | for(i=0;i<3;i++) this->values[3+i]=this->values[i];
|
---|
| 702 |
|
---|
| 703 | }
|
---|
| 704 | /*}}}*/
|
---|
[4471] | 705 | /*FUNCTION PentaVertexInput::VerticallyIntegrate{{{1*/
|
---|
| 706 | void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){
|
---|
| 707 |
|
---|
| 708 | /*Intermediaries*/
|
---|
| 709 | int i;
|
---|
| 710 | const int numgrids = 6;
|
---|
| 711 | int num_thickness_values;
|
---|
| 712 | double *thickness_values = NULL;
|
---|
| 713 |
|
---|
| 714 | /*Check that input provided is a thickness*/
|
---|
[5103] | 715 | if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumToString(thickness_input->EnumType()));
|
---|
[4471] | 716 |
|
---|
| 717 | /*Get Thickness value pointer*/
|
---|
| 718 | thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
|
---|
| 719 |
|
---|
| 720 | /*vertically integrate depending on type:*/
|
---|
| 721 | switch(thickness_input->Enum()){
|
---|
| 722 |
|
---|
| 723 | case PentaVertexInputEnum:
|
---|
| 724 | for(i=0;i<3;i++){
|
---|
| 725 | this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
|
---|
| 726 | this->values[i+3]=this->values[i];
|
---|
| 727 | }
|
---|
| 728 | return;
|
---|
| 729 |
|
---|
| 730 | default:
|
---|
| 731 | ISSMERROR("not implemented yet");
|
---|
| 732 | }
|
---|
| 733 | }
|
---|
| 734 | /*}}}*/
|
---|
| 735 | /*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
|
---|
| 736 | Input* PentaVertexInput::PointwiseDivide(Input* inputB){
|
---|
| 737 |
|
---|
| 738 | /*Ouput*/
|
---|
| 739 | PentaVertexInput* outinput=NULL;
|
---|
| 740 |
|
---|
| 741 | /*Intermediaries*/
|
---|
| 742 | int i;
|
---|
| 743 | PentaVertexInput *xinputB = NULL;
|
---|
| 744 | int B_numvalues;
|
---|
| 745 | double *B_values = NULL;
|
---|
| 746 | const int numgrids = 6;
|
---|
| 747 | double AdotBvalues[numgrids];
|
---|
| 748 |
|
---|
| 749 | /*Check that inputB is of the same type*/
|
---|
[5103] | 750 | if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumToString(inputB->Enum()));
|
---|
[4471] | 751 | xinputB=(PentaVertexInput*)inputB;
|
---|
| 752 |
|
---|
| 753 | /*Create point wise sum*/
|
---|
| 754 | for(i=0;i<numgrids;i++){
|
---|
| 755 | ISSMASSERT(xinputB->values[i]!=0);
|
---|
| 756 | AdotBvalues[i]=this->values[i]/xinputB->values[i];
|
---|
| 757 | }
|
---|
| 758 |
|
---|
[4899] | 759 | /*Create new Penta vertex input (copy of current input)*/
|
---|
[4471] | 760 | outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]);
|
---|
| 761 |
|
---|
| 762 | /*Return output pointer*/
|
---|
| 763 | return outinput;
|
---|
| 764 |
|
---|
| 765 | }
|
---|
| 766 | /*}}}*/
|
---|
[4546] | 767 | /*FUNCTION PentaVertexInput::GetVectorFromInputs{{{1*/
|
---|
[4048] | 768 | void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
|
---|
| 769 |
|
---|
| 770 | const int numvertices=6;
|
---|
[4502] | 771 | VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);
|
---|
[4048] | 772 |
|
---|
[4502] | 773 | } /*}}}*/
|
---|
[4546] | 774 | /*FUNCTION PentaVertexInput::GetValuesPtr{{{1*/
|
---|
[4057] | 775 | void PentaVertexInput::GetValuesPtr(double** pvalues,int* pnum_values){
|
---|
[4055] | 776 |
|
---|
| 777 | *pvalues=this->values;
|
---|
| 778 | *pnum_values=6;
|
---|
| 779 |
|
---|
| 780 | }
|
---|
| 781 | /*}}}*/
|
---|