Changeset 4832
- Timestamp:
- 07/27/10 11:41:32 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r4804 r4832 2125 2125 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/ 2126 2126 void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ 2127 2128 /*Inputs*/ 2129 bool collapsed,onbed; 2130 2131 /*Recover inputs*/ 2132 inputs->GetParameterValue(&collapsed,CollapseEnum); 2133 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2134 2135 /*MacAyeal, everything is done by the element on bed*/ 2136 if (collapsed){ 2137 if (!onbed){ 2138 return; 2139 } 2140 else{ 2141 InputUpdateFromSolutionDiagnosticMacAyeal(solution); 2142 return; 2143 } 2144 } 2145 else{ 2146 InputUpdateFromSolutionDiagnosticPattyn(solution); 2147 } 2148 } 2149 2150 /*}}}*/ 2151 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/ 2152 void Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){ 2153 2154 int i; 2155 2156 const int numvertices=6; 2157 const int numdofpervertex=2; 2158 const int numdof=numdofpervertex*numvertices; 2159 2160 int doflist[numdof]; 2161 double values[numdof]; 2162 double vx[numvertices]; 2163 double vy[numvertices]; 2164 double vz[numvertices]; 2165 double vel[numvertices]; 2166 int dummy; 2167 double pressure[numvertices]; 2168 double surface[numvertices]; 2169 double rho_ice,g; 2170 double xyz_list[numvertices][3]; 2171 double gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 2172 2173 Input *VzInput = NULL; 2174 double *VzPtr = NULL; 2175 Penta *penta = NULL; 2176 2177 /*OK, we are on bed. Now recover results*/ 2178 2179 /*Get dof list: */ 2180 GetDofList(&doflist[0],&dummy); 2181 2182 /*Get node data: */ 2183 GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); 2184 2185 /*Use the dof list to index into the solution vector: */ 2186 for(i=0;i<numdof;i++){ 2187 values[i]=solution[doflist[i]]; 2188 } 2189 2190 /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */ 2191 for(i=0;i<3;i++){ 2192 vx[i] =values[i*numdofpervertex+0]; 2193 vy[i] =values[i*numdofpervertex+1]; 2194 vx[i+3]=vx[i]; 2195 vy[i+3]=vy[i]; 2196 } 2197 2198 /*Get parameters fro pressure computation*/ 2199 rho_ice=matpar->GetRhoIce(); 2200 g=matpar->GetG(); 2201 2202 /*Start looping over all elements above current element and update all inputs*/ 2203 penta=this; 2204 for(;;){ 2205 2206 /*Now Compute vel*/ 2207 VzInput=inputs->GetInput(VzEnum); 2208 if (VzInput){ 2209 if (VzInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumAsString(VzInput->Enum())); 2210 VzInput->GetValuesPtr(&VzPtr,&dummy); 2211 for(i=0;i<numvertices;i++) vz[i]=VzPtr[i]; 2212 } 2213 else{for(i=0;i<numvertices;i++) vz[i]=0.0;} 2214 for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5); 2215 2216 /*Now compute pressure*/ 2217 inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); 2218 for(i=0;i<numvertices;i++){ 2219 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); 2220 } 2221 2222 /*Now, we have to move the previous Vx and Vy inputs to old 2223 * status, otherwise, we'll wipe them off: */ 2224 penta->inputs->ChangeEnum(VxEnum,VxOldEnum); 2225 penta->inputs->ChangeEnum(VyEnum,VyOldEnum); 2226 penta->inputs->ChangeEnum(PressureEnum,PressureOldEnum); 2227 2228 /*Add vx and vy as inputs to the tria element: */ 2229 penta->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); 2230 penta->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); 2231 penta->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 2232 penta->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 2233 2234 /*Stop if we have reached the surface*/ 2235 if (penta->IsOnSurface()) break; 2236 2237 /* get upper Penta*/ 2238 penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id); 2239 } 2240 } 2241 /*}}}*/ 2242 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{1*/ 2243 void Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){ 2127 2244 2128 2245 int i; … … 2145 2262 double gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 2146 2263 2147 Input * VzInput=NULL;2148 double * VzPtr=NULL;2264 Input *VzInput = NULL; 2265 double *VzPtr = NULL; 2149 2266 2150 2267 /*Get dof list: */ … … 2186 2303 g=matpar->GetG(); 2187 2304 inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); 2188 2305 2189 2306 for(i=0;i<numvertices;i++){ 2190 2307 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); … … 2199 2316 this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); 2200 2317 this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); 2201 this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));2318 this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 2202 2319 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 2203 2320 } … … 4916 5033 4917 5034 /*Are we on the base, not on the surface?:*/ 4918 if(onbed ==1){5035 if(onbed){ 4919 5036 4920 5037 /*OK, we are on bed. we will follow the steps: -
issm/trunk/src/c/objects/Elements/Penta.h
r4780 r4832 177 177 void InputUpdateFromSolutionBalancedvelocities( double* solutiong); 178 178 void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); 179 void InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong); 180 void InputUpdateFromSolutionDiagnosticPattyn( double* solutiong); 179 181 void InputUpdateFromSolutionDiagnosticHutter( double* solutiong); 180 182 void InputUpdateFromSolutionDiagnosticVert( double* solutiong); -
issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
r4592 r4832 128 128 } 129 129 130 131 /*extrude if we are in 3D: */132 if (dim==3){133 if(verbose)_printf_("%s\n","extruding velocity of collapsed elements...");134 InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,true); //true means that the velocity is only extruded if collapsed135 InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,true);136 }137 138 130 //more output might be needed, when running in control 139 131 if(pKff0){
Note:
See TracChangeset
for help on using the changeset viewer.