![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 16947
- Timestamp:
- 11/26/13 08:36:55 (11 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r16936 r16947 83 83 } 84 84 /*}}}*/ 85 void Element::GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity){/*{{{*/ 86 /*Compute deformational heating from epsilon and viscosity */ 87 88 IssmDouble epsilon_matrix[3][3]; 89 IssmDouble epsilon_eff; 90 IssmDouble epsilon_sqr[3][3]; 91 92 /* Build epsilon matrix */ 93 epsilon_matrix[0][0]=epsilon[0]; 94 epsilon_matrix[1][0]=epsilon[3]; 95 epsilon_matrix[2][0]=epsilon[4]; 96 epsilon_matrix[0][1]=epsilon[3]; 97 epsilon_matrix[1][1]=epsilon[1]; 98 epsilon_matrix[2][1]=epsilon[5]; 99 epsilon_matrix[0][2]=epsilon[4]; 100 epsilon_matrix[1][2]=epsilon[5]; 101 epsilon_matrix[2][2]=epsilon[2]; 102 103 /* Effective value of epsilon_matrix */ 104 epsilon_sqr[0][0]=epsilon_matrix[0][0]*epsilon_matrix[0][0]; 105 epsilon_sqr[1][0]=epsilon_matrix[1][0]*epsilon_matrix[1][0]; 106 epsilon_sqr[2][0]=epsilon_matrix[2][0]*epsilon_matrix[2][0]; 107 epsilon_sqr[0][1]=epsilon_matrix[0][1]*epsilon_matrix[0][1]; 108 epsilon_sqr[1][1]=epsilon_matrix[1][1]*epsilon_matrix[1][1]; 109 epsilon_sqr[2][1]=epsilon_matrix[2][1]*epsilon_matrix[2][1]; 110 epsilon_sqr[0][2]=epsilon_matrix[0][2]*epsilon_matrix[0][2]; 111 epsilon_sqr[1][2]=epsilon_matrix[1][2]*epsilon_matrix[1][2]; 112 epsilon_sqr[2][2]=epsilon_matrix[2][2]*epsilon_matrix[2][2]; 113 epsilon_eff=1/sqrt(2.)*sqrt(epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]); 114 115 /*Phi = Tr(sigma * eps) 116 * = Tr(sigma'* eps) 117 * = 2 * eps_eff * sigma'_eff 118 * = 4 * mu * eps_eff ^2*/ 119 *phi=4.*epsilon_eff*epsilon_eff*viscosity; 120 } 121 /*}}}*/ 85 122 Input* Element::GetInput(int inputenum){/*{{{*/ 86 123 return inputs->GetInput(inputenum); 87 124 }/*}}}*/ 125 void Element::GetVerticesSidList(int* sidlist){/*{{{*/ 126 127 int numvertices = this->GetNumberOfVertices(); 128 for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid(); 129 } 130 /*}}}*/ 131 void Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/ 132 133 int numvertices = this->GetNumberOfVertices(); 134 for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity(); 135 } 136 /*}}}*/ 88 137 bool Element::IsFloating(){/*{{{*/ 89 138 … … 104 153 return shelf; 105 154 }/*}}}*/ 155 void Element::ResultInterpolation(int* pinterpolation,int output_enum){/*{{{*/ 156 157 Input* input=this->inputs->GetInput(output_enum); 158 159 /*If this input is not already in Inputs, maybe it needs to be computed?*/ 160 if(!input){ 161 switch(output_enum){ 162 case ViscousHeatingEnum: 163 this->ViscousHeatingCreateInput(); 164 input=this->inputs->GetInput(output_enum); 165 break; 166 case StressTensorxxEnum: 167 this->ComputeStressTensor(); 168 input=this->inputs->GetInput(output_enum); 169 break; 170 case StressTensorxyEnum: 171 this->ComputeStressTensor(); 172 input=this->inputs->GetInput(output_enum); 173 break; 174 case StressTensorxzEnum: 175 this->ComputeStressTensor(); 176 input=this->inputs->GetInput(output_enum); 177 break; 178 case StressTensoryyEnum: 179 this->ComputeStressTensor(); 180 input=this->inputs->GetInput(output_enum); 181 break; 182 case StressTensoryzEnum: 183 this->ComputeStressTensor(); 184 input=this->inputs->GetInput(output_enum); 185 break; 186 case StressTensorzzEnum: 187 this->ComputeStressTensor(); 188 input=this->inputs->GetInput(output_enum); 189 break; 190 default: 191 _error_("input "<<EnumToStringx(output_enum)<<" not found in element"); 192 } 193 } 194 195 /*Assign output pointer*/ 196 *pinterpolation = input->GetResultInterpolation(); 197 }/*}}}*/ 198 void Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/ 199 200 Input* input=this->inputs->GetInput(output_enum); 201 if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element"); 202 203 switch(input->GetResultInterpolation()){ 204 case P0Enum:{ 205 IssmDouble value; 206 bool bvalue; 207 Input* input = this->GetInput(output_enum); _assert_(input); 208 switch(input->ObjectEnum()){ 209 case DoubleInputEnum: 210 input->GetInputValue(&value); 211 break; 212 case BoolInputEnum: 213 input->GetInputValue(&bvalue); 214 value=reCast<IssmDouble>(bvalue); 215 break; 216 default: 217 _error_("Input of type "<<EnumToStringx(input->GetResultInterpolation())<<" not supported"); 218 } 219 vector->SetValue(this->Sid(),value,INS_VAL); 220 break; 221 } 222 case P1Enum:{ 223 int numvertices = this->GetNumberOfVertices(); 224 IssmDouble *values = xNew<IssmDouble>(numvertices); 225 int *connectivity= xNew<int>(numvertices); 226 int *sidlist = xNew<int>(numvertices); 227 228 this->GetVerticesSidList(sidlist); 229 this->GetVerticesConnectivityList(connectivity); 230 this->GetInputListOnVertices(values,output_enum); 231 for(int i=0;i<numvertices;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]); 232 233 vector->SetValues(numvertices,sidlist,values,ADD_VAL); 234 235 xDelete<IssmDouble>(values); 236 xDelete<int>(connectivity); 237 xDelete<int>(sidlist); 238 break; 239 } 240 default: 241 _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet"); 242 } 243 } /*}}}*/ 106 244 IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/ 107 245 _assert_(matpar); 108 246 return this->matpar->TMeltingPoint(pressure); 109 247 }/*}}}*/ 110 248 void Element::ViscousHeatingCreateInput(void){/*{{{*/ 249 250 /*Intermediaries*/ 251 IssmDouble phi; 252 IssmDouble viscosity; 253 IssmDouble epsilon[6]; 254 IssmDouble thickness; 255 IssmDouble *xyz_list = NULL; 256 257 /*Fetch number vertices and allocate memory*/ 258 int numvertices = this->GetNumberOfVertices(); 259 IssmDouble* viscousheating = xNew<IssmDouble>(numvertices); 260 261 /*Retrieve all inputs and parameters*/ 262 this->GetVerticesCoordinatesBase(&xyz_list); 263 Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); 264 Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); 265 Input* vz_input = this->GetInput(VzEnum); _assert_(vz_input); 266 Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input); 267 268 /*loop over vertices: */ 269 Gauss* gauss=this->NewGauss(); 270 for (int iv=0;iv<numvertices;iv++){ 271 gauss->GaussVertex(iv); 272 273 thickness_input->GetInputValue(&thickness,gauss); 274 275 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); 276 this->material->GetViscosity3dFS(&viscosity,&epsilon[0]); 277 this->GetPhi(&phi,&epsilon[0],viscosity); 278 279 viscousheating[iv]=phi*thickness; 280 } 281 282 /*Create PentaVertex input, which will hold the basal friction:*/ 283 this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum); 284 285 /*Clean up and return*/ 286 xDelete<IssmDouble>(viscousheating); 287 delete gauss; 288 } 289 /*}}}*/ 111 290 void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ 112 291 -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r16936 r16947 61 61 Input* GetInput(int inputenum); 62 62 IssmDouble GetMaterialParameter(int enum_in); 63 void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); 64 void GetVerticesSidList(int* sidlist); 65 void GetVerticesConnectivityList(int* connectivitylist); 63 66 bool IsFloating(); 67 void ResultInterpolation(int* pinterpolation,int output_enum); 68 void ResultToVector(Vector<IssmDouble>* vector,int output_enum); 64 69 void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); 65 70 void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); 66 71 void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); 67 72 IssmDouble TMeltingPoint(IssmDouble pressure); 73 void ViscousHeatingCreateInput(void); 68 74 void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); 69 75 void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); … … 77 83 void TransformLoadVectorCoord(ElementVector* pe,int transformenum); 78 84 void TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list); 79 void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");}; 80 void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");}; 85 void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};/*Tiling only*/ 86 void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/ 81 87 void TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* transformenum_list); 82 88 void TransformSolutionCoord(IssmDouble* values,int transformenum); … … 87 93 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list); 88 94 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum); 89 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");}; 95 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/ 90 96 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list); 91 97 … … 168 174 virtual void ComputeStrainRate(Vector<IssmDouble>* eps)=0; 169 175 virtual void ComputeStressTensor(void)=0; 170 virtual void ResultInterpolation(int* pinterpolation,int output_enum)=0; 171 virtual void ResultToVector(Vector<IssmDouble>* vector,int output_enum)=0; 176 172 177 virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0; 173 178 virtual void InputDuplicate(int original_enum,int new_enum)=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r16936 r16947 1153 1153 } 1154 1154 /*}}}*/ 1155 /*FUNCTION Penta::GetVertexSidList{{{*/1156 void Penta::GetVertexSidList(int* sidlist){1157 1158 for(int i=0;i<NUMVERTICES;i++) sidlist[i]=vertices[i]->Sid();1159 1160 }1161 /*}}}*/1162 /*FUNCTION Penta::GetConnectivityList {{{*/1163 void Penta::GetConnectivityList(int* connectivity){1164 for(int i=0;i<NUMVERTICES;i++) connectivity[i]=vertices[i]->Connectivity();1165 }1166 /*}}}*/1167 1155 /*FUNCTION Penta::GetElementType {{{*/ 1168 1156 int Penta::GetElementType(){ … … 1448 1436 input->GetInputValue(pvalue,gauss); 1449 1437 delete gauss; 1450 }1451 /*}}}*/1452 /*FUNCTION Penta::GetPhi {{{*/1453 void Penta::GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity){1454 /*Compute deformational heating from epsilon and viscosity */1455 1456 IssmDouble epsilon_matrix[3][3];1457 IssmDouble epsilon_eff;1458 IssmDouble epsilon_sqr[3][3];1459 1460 /* Build epsilon matrix */1461 epsilon_matrix[0][0]=*(epsilon+0);1462 epsilon_matrix[1][0]=*(epsilon+3);1463 epsilon_matrix[2][0]=*(epsilon+4);1464 epsilon_matrix[0][1]=*(epsilon+3);1465 epsilon_matrix[1][1]=*(epsilon+1);1466 epsilon_matrix[2][1]=*(epsilon+5);1467 epsilon_matrix[0][2]=*(epsilon+4);1468 epsilon_matrix[1][2]=*(epsilon+5);1469 epsilon_matrix[2][2]=*(epsilon+2);1470 1471 /* Effective value of epsilon_matrix */1472 epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2);1473 epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2);1474 epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2);1475 epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2);1476 epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2);1477 epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2);1478 epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2);1479 epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2);1480 epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);1481 epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5);1482 1483 /*Phi = Tr(sigma * eps)1484 * = Tr(sigma'* eps)1485 * = 2 * eps_eff * sigma'_eff1486 * = 4 * mu * eps_eff ^2*/1487 *phi=4*pow(epsilon_eff,2.0)*viscosity;1488 1438 } 1489 1439 /*}}}*/ … … 2724 2674 } 2725 2675 } 2726 }2727 /*}}}*/2728 /*FUNCTION Penta::ResultInterpolation{{{*/2729 void Penta::ResultInterpolation(int* pinterpolation,int output_enum){2730 2731 Input* input=this->inputs->GetInput(output_enum);2732 2733 /*If this input is not already in Inputs, maybe it needs to be computed?*/2734 if(!input){2735 switch(output_enum){2736 case ViscousHeatingEnum:2737 this->ViscousHeatingCreateInput();2738 input=this->inputs->GetInput(output_enum);2739 break;2740 case StressTensorxxEnum:2741 this->ComputeStressTensor();2742 input=this->inputs->GetInput(output_enum);2743 break;2744 case StressTensorxyEnum:2745 this->ComputeStressTensor();2746 input=this->inputs->GetInput(output_enum);2747 break;2748 case StressTensorxzEnum:2749 this->ComputeStressTensor();2750 input=this->inputs->GetInput(output_enum);2751 break;2752 case StressTensoryyEnum:2753 this->ComputeStressTensor();2754 input=this->inputs->GetInput(output_enum);2755 break;2756 case StressTensoryzEnum:2757 this->ComputeStressTensor();2758 input=this->inputs->GetInput(output_enum);2759 break;2760 case StressTensorzzEnum:2761 this->ComputeStressTensor();2762 input=this->inputs->GetInput(output_enum);2763 break;2764 default:2765 _error_("input "<<EnumToStringx(output_enum)<<" not found in element");2766 }2767 }2768 2769 /*Assign output pointer*/2770 *pinterpolation = input->GetResultInterpolation();2771 2772 }2773 /*}}}*/2774 /*FUNCTION Penta::ResultToVector{{{*/2775 void Penta::ResultToVector(Vector<IssmDouble>* vector,int output_enum){2776 2777 Input* input=this->inputs->GetInput(output_enum);2778 if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");2779 2780 switch(input->GetResultInterpolation()){2781 case P0Enum:2782 _error_("P0 not implemented yet for input "<<EnumToStringx(output_enum));2783 break;2784 case P1Enum:{2785 IssmDouble values[NUMVERTICES];2786 int connectivity[NUMVERTICES];2787 int sidlist[NUMVERTICES];2788 2789 this->GetVertexSidList(&sidlist[0]);2790 this->GetConnectivityList(&connectivity[0]);2791 this->GetInputListOnVertices(&values[0],output_enum);2792 for(int i=0;i<NUMVERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);2793 2794 vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);2795 break;2796 }2797 default:2798 _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");2799 }2800 2801 2802 2676 } 2803 2677 /*}}}*/ … … 3356 3230 } 3357 3231 /*}}}*/ 3358 /*FUNCTION Penta::ViscousHeatingCreateInput {{{*/3359 void Penta::ViscousHeatingCreateInput(void){3360 3361 /*Intermediaries*/3362 IssmDouble phi;3363 IssmDouble viscosity;3364 IssmDouble xyz_list[NUMVERTICES][3];3365 IssmDouble epsilon[6];3366 IssmDouble viscousheating[NUMVERTICES];3367 IssmDouble thickness;3368 3369 /*Retrieve all inputs and parameters*/3370 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3371 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3372 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3373 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);3374 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);3375 3376 /*loop over vertices: */3377 GaussPenta* gauss=new GaussPenta();3378 for (int iv=0;iv<NUMVERTICES;iv++){3379 gauss->GaussVertex(iv);3380 3381 thickness_input->GetInputValue(&thickness,gauss);3382 3383 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);3384 material->GetViscosity3dFS(&viscosity,&epsilon[0]);3385 GetPhi(&phi, &epsilon[0], viscosity);3386 3387 viscousheating[iv]=phi*thickness;3388 }3389 3390 /*Create PentaVertex input, which will hold the basal friction:*/3391 this->inputs->AddInput(new PentaInput(ViscousHeatingEnum,&viscousheating[0],P1Enum));3392 3393 /*Clean up and return*/3394 delete gauss;3395 }3396 /*}}}*/3397 3232 /*FUNCTION Penta::ViscousHeating{{{*/ 3398 3233 void Penta::ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ … … 5058 4893 IssmDouble value,gradient; 5059 4894 5060 this->Get ConnectivityList(&connectivity[0]);5061 this->GetVert exSidList(&sidlist[0]);4895 this->GetVerticesConnectivityList(&connectivity[0]); 4896 this->GetVerticesSidList(&sidlist[0]); 5062 4897 5063 4898 GaussPenta* gauss=new GaussPenta(); -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r16926 r16947 108 108 bool IsZeroLevelset(int levelset_enum); 109 109 void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); 110 111 void ResultInterpolation(int* pinterpolation,int output_enum);112 void ResultToVector(Vector<IssmDouble>* vector,int output_enum);113 110 void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); 114 111 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); … … 125 122 int VertexConnectivity(int vertexindex); 126 123 void VerticalSegmentIndices(int** pindices,int* pnumseg); 127 void ViscousHeatingCreateInput(void);128 124 void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); 129 125 … … 213 209 214 210 void GetVertexPidList(int* doflist); 215 void GetVertexSidList(int* sidlist);216 void GetConnectivityList(int* connectivity);217 211 int GetElementType(void); 218 212 Input* GetMaterialInput(int inputenum); … … 228 222 void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); 229 223 Node* GetNode(int node_number); 230 void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity);231 224 void GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input); 232 225 void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input); -
issm/trunk-jpl/src/c/classes/Elements/Seg.h
r16926 r16947 179 179 int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){_error_("not implemented yet");}; 180 180 void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; 181 void ResultInterpolation(int* pinterpolation,int output_enum){_error_("not implemented");};182 void ResultToVector(Vector<IssmDouble>* vector,int output_enum){_error_("not implemented");};183 181 void ResetCoordinateSystem(void){_error_("not implemented yet");}; 184 182 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r16916 r16947 1510 1510 } 1511 1511 /*}}}*/ 1512 /*FUNCTION Tria::GetVertexSidList {{{*/1513 void Tria::GetVertexSidList(int* sidlist){1514 for(int i=0;i<NUMVERTICES;i++) sidlist[i]=vertices[i]->Sid();1515 }1516 /*}}}*/1517 /*FUNCTION Tria::GetConnectivityList {{{*/1518 void Tria::GetConnectivityList(int* connectivity){1519 for(int i=0;i<NUMVERTICES;i++) connectivity[i]=vertices[i]->Connectivity();1520 }1521 /*}}}*/1522 1512 /*FUNCTION Tria::GetVectorFromInputs{{{*/ 1523 1513 void Tria::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){ … … 2537 2527 } 2538 2528 /*}}}*/ 2539 /*FUNCTION Tria::ResultInterpolation{{{*/2540 void Tria::ResultInterpolation(int* pinterpolation,int output_enum){2541 2542 Input* input=this->inputs->GetInput(output_enum);2543 2544 /*If this input is not already in Inputs, maybe it needs to be computed?*/2545 if(!input){2546 switch(output_enum){2547 case StressTensorxxEnum:2548 this->ComputeStressTensor();2549 input=this->inputs->GetInput(output_enum);2550 break;2551 case StressTensorxyEnum:2552 this->ComputeStressTensor();2553 input=this->inputs->GetInput(output_enum);2554 break;2555 case StressTensorxzEnum:2556 this->ComputeStressTensor();2557 input=this->inputs->GetInput(output_enum);2558 break;2559 case StressTensoryyEnum:2560 this->ComputeStressTensor();2561 input=this->inputs->GetInput(output_enum);2562 break;2563 case StressTensoryzEnum:2564 this->ComputeStressTensor();2565 input=this->inputs->GetInput(output_enum);2566 break;2567 case StressTensorzzEnum:2568 this->ComputeStressTensor();2569 input=this->inputs->GetInput(output_enum);2570 break;2571 default:2572 _error_("input "<<EnumToStringx(output_enum)<<" not found in element");2573 }2574 }2575 2576 /*Assign output pointer*/2577 *pinterpolation = input->GetResultInterpolation();2578 2579 }2580 /*}}}*/2581 /*FUNCTION Tria::ResultToVector{{{*/2582 void Tria::ResultToVector(Vector<IssmDouble>* vector,int output_enum){2583 2584 Input* input=this->inputs->GetInput(output_enum);2585 if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");2586 2587 switch(input->GetResultInterpolation()){2588 case P0Enum:2589 _error_("not implemented...");2590 break;2591 case P1Enum:{2592 IssmDouble values[NUMVERTICES];2593 int connectivity[NUMVERTICES];2594 int sidlist[NUMVERTICES];2595 2596 this->GetVertexSidList(&sidlist[0]);2597 this->GetConnectivityList(&connectivity[0]);2598 this->GetInputListOnVertices(&values[0],output_enum);2599 for(int i=0;i<NUMVERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);2600 2601 vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);2602 break;2603 }2604 default:2605 _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");2606 }2607 2608 2609 }2610 /*}}}*/2611 2529 /*FUNCTION Tria::SetClone {{{*/ 2612 2530 void Tria::SetClone(int* minranks){ … … 2997 2915 2998 2916 /*Figure out the average for this element: */ 2999 this->GetVert exSidList(&offsetsid[0]);2917 this->GetVerticesSidList(&offsetsid[0]); 3000 2918 this->GetVertexPidList(&offsetdof[0]); 3001 2919 mean=0; … … 4569 4487 IssmDouble value,gradient; 4570 4488 4571 this->Get ConnectivityList(&connectivity[0]);4572 this->GetVert exSidList(&sidlist[0]);4489 this->GetVerticesConnectivityList(&connectivity[0]); 4490 this->GetVerticesSidList(&sidlist[0]); 4573 4491 4574 4492 GaussTria* gauss=new GaussTria(); … … 4832 4750 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 4833 4751 GradientIndexing(&vertexpidlist[0],control_index); 4834 this->Get ConnectivityList(&connectivity[0]);4752 this->GetVerticesConnectivityList(&connectivity[0]); 4835 4753 4836 4754 /*Build frictoin element, needed later: */ -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r16926 r16947 113 113 int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); 114 114 void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); 115 void ResultInterpolation(int* pinterpolation,int output_enum);116 void ResultToVector(Vector<IssmDouble>* vector,int output_enum);117 115 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); 118 116 void ResetCoordinateSystem(void); … … 247 245 248 246 void GetVertexPidList(int* doflist); 249 void GetVertexSidList(int* sidlist);250 void GetConnectivityList(int* connectivity);251 247 IssmDouble GetYcoord(Gauss* gauss); 252 248 IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented");};
Note:
See TracChangeset
for help on using the changeset viewer.
