 |
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
|
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
00001 /************************************************************************ 00002 * MechSys - Open Library for Mechanical Systems * 00003 * Copyright (C) 2005 Dorival M. Pedroso * 00004 * * 00005 * This program is free software: you can redistribute it and/or modify * 00006 * it under the terms of the GNU General Public License as published by * 00007 * the Free Software Foundation, either version 3 of the License, or * 00008 * any later version. * 00009 * * 00010 * This program is distributed in the hope that it will be useful, * 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00013 * GNU General Public License for more details. * 00014 * * 00015 * You should have received a copy of the GNU General Public License * 00016 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00017 ************************************************************************/ 00018 00019 #ifndef MPM_PROBLEMS2D_H 00020 #define MPM_PROBLEMS2D_H 00021 00022 // FLTK 00023 #include <FL/Fl_Choice.H> 00024 00025 // MechSys 00026 #include <mechsys/util/array.h> 00027 #include<mechsys/models/linelastic.h> 00028 #include<mechsys/models/neohookean.h> 00029 #include<mechsys/models/camclay.h> 00030 #include<mechsys/models/elastoplastic.h> 00031 00032 // Local 00033 #include <mechsys/mpm/defs.h> 00034 #include <mechsys/mpm/grid2d.h> 00035 #include <mechsys/mpm/mpoints2d.h> 00036 00037 namespace MPM { 00038 00039 // Fixed nodes? 00040 bool CbIsFixed1 (Vector3D const & N, FixType & FType) 00041 { 00042 if (N(0)<1.01) { FType = FIX_XY; return true; } 00043 if (N(1)<0.11 || N(1)>0.19) { FType = FIX_Y; return true; } 00044 return false; 00045 } 00046 bool CbIsFixed2 (Vector3D const & N, FixType & FType) 00047 { 00048 if (N(0)<0.01 && N(1)<0.01) { FType=FIX_XY; return true; } 00049 if (N(0)<0.01) { FType=FIX_X; return true; } 00050 if (N(1)<0.01) { FType=FIX_Y; return true; } 00051 return false; 00052 } 00053 bool CbIsFixed3 (Vector3D const & N, FixType & FType) 00054 { 00055 /* 00056 if (N(0)>-0.051 && N(0)<-0.049 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_XY; return true; } 00057 if (N(0)> 0.99 && N(0)< 1.01 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_XY; return true; } 00058 if (N(0)> 0.99 && N(0)< 1.01 && N(1)> 0.99 && N(1)<1.01) { FType=FIX_XY; return true; } 00059 if (N(0)>-0.051 && N(0)<-0.049 && N(1)> 0.99 && N(1)<1.01) { FType=FIX_XY; return true; } 00060 if (N(0)>-0.051 && N(0)<-0.049 && N(1)>-0.049 && N(1)<1.01) { FType=FIX_X; return true; } 00061 if (N(0)< 1.01 && N(0)>-0.051 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_Y; return true; } 00062 if (N(0)> 0.99 && N(0)< 1.01 && N(1)>-0.049 && N(1)<1.01) { FType=FIX_X; return true; } 00063 if (N(0)< 1.01 && N(0)>-0.051 && N(1)> 0.99 && N(1)<1.01) { FType=FIX_Y; return true; } 00064 */ 00065 if (N(0)>-0.051 && N(0)<-0.049 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_XY; return true; } 00066 if (N(0)> 0.69 && N(0)< 0.71 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_XY; return true; } 00067 if (N(0)> 0.69 && N(0)< 0.71 && N(1)> 0.69 && N(1)<0.71) { FType=FIX_XY; return true; } 00068 if (N(0)>-0.051 && N(0)<-0.049 && N(1)> 0.69 && N(1)<0.71) { FType=FIX_XY; return true; } 00069 if (N(0)>-0.051 && N(0)<-0.049 && N(1)>-0.049 && N(1)<0.71) { FType=FIX_X; return true; } 00070 if (N(0)< 0.71 && N(0)>-0.051 && N(1)>-0.051 && N(1)<-0.049) { FType=FIX_Y; return true; } 00071 if (N(0)> 0.69 && N(0)< 0.71 && N(1)>-0.049 && N(1)<0.71) { FType=FIX_X; return true; } 00072 if (N(0)< 0.71 && N(0)>-0.051 && N(1)> 0.69 && N(1)<0.71) { FType=FIX_Y; return true; } 00073 return false; 00074 } 00075 bool CbIsFixed4 (Vector3D const & N, FixType & FType) 00076 { 00077 if (N(0)<0.01 && N(1)<0.01) { FType=FIX_XY; return true; } 00078 if (N(0)<0.01) { FType=FIX_X; return true; } 00079 if (N(1)<0.01) { FType=FIX_Y; return true; } 00080 return false; 00081 } 00082 bool CbIsFixed5 (Vector3D const & N, FixType & FType) // micrometer 00083 { 00084 if (N(0)>-0.1 && N(0)<0.1) { FType=FIX_X; return true; } 00085 return false; 00086 } 00087 bool CbIsFixed6 (Vector3D const & N, FixType & FType) 00088 { 00089 /* 00090 if (N(0)>-0.01 && N(0)<0.01 && N(1)>-0.01 && N(1)<0.01) { FType=FIX_XY; return true; } 00091 if (N(0)> 19.99 && N(0)< 20.01 && N(1)>-0.01 && N(1)<0.01) { FType=FIX_XY; return true; } 00092 if (N(0)> 19.99 && N(0)< 20.01 && N(1)> 19.99 && N(1)<20.01) { FType=FIX_XY; return true; } 00093 if (N(0)>-0.01 && N(0)<0.01 && N(1)> 19.99 && N(1)<20.01) { FType=FIX_XY; return true; } 00094 if (N(0)>-0.01 && N(0)<0.01 && N(1)>0.01 && N(1)<20.01) { FType=FIX_X; return true; } 00095 if (N(0)< 20.01 && N(0)>-0.01 && N(1)>-0.01 && N(1)<0.01) { FType=FIX_Y; return true; } 00096 if (N(0)> 19.99 && N(0)< 20.01 && N(1)>0.01 && N(1)<20.01) { FType=FIX_X; return true; } 00097 if (N(0)< 20.01 && N(0)>-0.01 && N(1)> 19.99 && N(1)<20.01) { FType=FIX_Y; return true; } 00098 */ 00099 return false; 00100 } 00101 bool CbIsFixed7 (Vector3D const & N, FixType & FType) 00102 { 00103 if (N(0)<0.01) { FType=FIX_XY; return true; } 00104 else return false; 00105 } 00106 bool CbIsFixed8 (Vector3D const & N, FixType & FType) 00107 { 00108 if (N(1)<0.01) { FType=FIX_XY; return true; } 00109 else return false; 00110 } 00111 bool CbIsFixed10(Vector3D const & N, FixType & FType) 00112 { 00113 if (N(0)<0.01) { FType = FIX_XY; return true; } 00114 else { FType = FIX_Y; return true; } 00115 } 00116 bool CbIsFixed11(Vector3D const & N, FixType & FType) 00117 { 00118 if (N(0)<0.01 && N(1)<0.01) { FType = FIX_XY; return true; } 00119 else if (N(0)<0.01) { FType = FIX_X; return true; } 00120 else if (N(1)<0.01) { FType = FIX_Y; return true; } 00121 return false; 00122 } 00123 bool CbIsFixed12(Vector3D const & N, FixType & FType) 00124 { 00125 if (N(1)<0.01) // bottom 00126 { 00127 if (fabs(N(0)-20.0)<0.01) { FType = FIX_XY; return true; } 00128 else { FType = FIX_Y; return true; } 00129 } 00130 return false; 00131 } 00132 bool CbIsFixed13(Vector3D const & N, FixType & FType) 00133 { 00134 if (N(1)<0.01) // bottom 00135 { 00136 if (N(0)<0.01 || N(0)>14.95) { FType = FIX_XY; return true; } 00137 else { FType = FIX_Y; return true; } 00138 } 00139 else if (N(0)<0.01 || N(0)>14.95) // left or right 00140 { 00141 FType = FIX_X; return true; 00142 } 00143 return false; 00144 } 00145 bool CbIsFixed14(Vector3D const & N, FixType & FType) 00146 { 00147 double L = 30.0; 00148 if (N(1)<0.01) // bottom 00149 { 00150 if (N(0)<0.01 || N(0)>L-0.5) { FType = FIX_XY; return true; } 00151 else { FType = FIX_Y; return true; } 00152 } 00153 else if (N(0)<0.01 || N(0)>L-0.5) // left or right 00154 { 00155 FType = FIX_X; return true; 00156 } 00157 return false; 00158 } 00159 bool CbIsFixed15(Vector3D const & N, FixType & FType) 00160 { 00161 if (fabs(N(0)-5.0)<0.01 && fabs(N(1)-5.0)<0.01) 00162 { 00163 FType = FIX_XY; 00164 return true; 00165 } 00166 return false; 00167 } 00168 bool CbIsFixed16(Vector3D const & N, FixType & FType) 00169 { 00170 if (fabs(N(1)-9.5)<0.05) 00171 { 00172 FType = FIX_Y; 00173 return true; 00174 } 00175 return false; 00176 } 00177 bool CbIsFixed17(Vector3D const & N, FixType & FType) 00178 { 00179 if (N(1)>9.4) 00180 { 00181 if (fabs(N(0)-1.0)<0.01) FType = FIX_XY; 00182 else FType = FIX_Y; 00183 return true; 00184 } 00185 return false; 00186 } 00187 00188 // Is point inside geometry callbacks 00189 bool CbIsPointInGeom1 (Vector3D const & P, int & ClrIdx) 00190 { 00191 00192 if (P(0)>=1.0 && P(0)<=2.0 && P(1)>=1.50 && P(1)<=2.0) return true; 00193 else return false; 00194 00195 } 00196 bool CbIsPointInGeom2 (Vector3D const & P, int & ClrIdx) 00197 { 00198 if (P(0)>0.0 && P(0)<1.0 && P(1)>0.0 && P(1)<1.0) return true; 00199 else return false; 00200 } 00201 bool CbIsPointInGeom3 (Vector3D const & P, int & ClrIdx) 00202 { 00203 // Cylinders 00204 /* 00205 double xc1=0.25; double yc1=0.25; double r1=0.14; 00206 double xc2=0.70; double yc2=0.70; double r2=0.14; 00207 */ 00208 double xc1=0.15; double yc1=0.15; double r1=0.14; 00209 double xc2=0.50; double yc2=0.50; double r2=0.14; 00210 00211 // Bottom-left cylinder 00212 double d1 = sqrt(pow(P(0)-xc1,2.0)+pow(P(1)-yc1,2.0)); 00213 00214 // Upper-right cylinder 00215 double d2 = sqrt(pow(P(0)-xc2,2.0)+pow(P(1)-yc2,2.0)); 00216 00217 if (d1<=r1 || d2<=r2) return true; 00218 else return false; 00219 /* 00220 else 00221 { 00222 if (P(0)<0.05 && P(0)>0 && P(1)>0 && P(1)<0.95) return true; // left 00223 if (P(0)>0.9 && P(0)<0.95 && P(1)>0 && P(1)<0.95) return true; // right 00224 if (P(0)>0 && P(0)<0.95 && P(1)<0.05 && P(1)>0 ) return true; // bottom 00225 if (P(0)>0 && P(0)<0.95 && P(1)>0.9 && P(1)<0.95) return true; // top 00226 return false; 00227 } 00228 */ 00229 } 00230 bool CbIsPointInGeom4 (Vector3D const & P, int & ClrIdx) 00231 { 00232 if (P(0)>0.0 && P(0)<1.0 && P(1)>0.0 && P(1)<1.0) return true; 00233 else return false; 00234 } 00235 bool CbIsPointInGeom5 (Vector3D const & P, int & ClrIdx) // micrometer 00236 { 00237 if (P(0)>0.0 && P(0)<60.0 && P(1)>0.0 && P(1)<40.0) return true; 00238 else return false; 00239 } 00240 bool CbIsPointInGeom6 (Vector3D const & P, int & ClrIdx) 00241 { 00242 // Cylinders 00243 double xc=10.0; double yc=10.0; double r=2.5; 00244 00245 // Distance 00246 double d = sqrt(pow(P(0)-xc,2.0)+pow(P(1)-yc,2.0)); 00247 00248 if (d<=r) return true; 00249 else 00250 { 00251 if (P(0)> 2.0 && P(0)< 3.0 && P(1)> 2.0 && P(1)<18.0) return true; 00252 if (P(0)>17.0 && P(0)<18.0 && P(1)> 2.0 && P(1)<18.0) return true; 00253 if (P(0)> 2.0 && P(0)<18.0 && P(1)>17.0 && P(1)<18.0) return true; 00254 if (P(0)> 2.0 && P(0)<18.0 && P(1)> 2.0 && P(1)< 3.0) return true; 00255 return false; 00256 } 00257 } 00258 bool CbIsPointInGeom7 (Vector3D const & P, int & ClrIdx) 00259 { 00260 if (P(0)>0.0 && P(0)<19.0 && P(1)>10.0 && P(1)<12.0) return true; 00261 else return false; 00262 } 00263 bool CbIsPointInGeom8 (Vector3D const & P, int & ClrIdx) 00264 { 00265 if (P(0)>5.0 && P(0)<7.0 && P(1)>0.0 && P(1)<12.0) return true; 00266 else return false; 00267 } 00268 bool CbIsPointInGeom10(Vector3D const & P, int & ClrIdx) 00269 { 00270 if (P(0)>=0.0 && P(0)<=25.0 && P(1)>=1.0 && P(1)<=2.0) return true; 00271 else return false; 00272 } 00273 bool CbIsPointInGeom11(Vector3D const & P, int & ClrIdx) 00274 { 00275 if (P(0)>0.0 && P(0)<10.01 && P(1)>0.0 && P(1)<10.01) 00276 { 00277 double xc1=0.0; double yc1=0.0; double r1=1.0; 00278 double d1 = sqrt(pow(P(0)-xc1,2.0)+pow(P(1)-yc1,2.0)); 00279 if (d1>=r1) return true; 00280 } 00281 return false; 00282 } 00283 bool CbIsPointInGeom12(Vector3D const & P, int & ClrIdx) 00284 { 00285 if (P(1)<=30.0-P(0) && P(1)<=P(0)-10.0 && P(1)>=0.0) return true; 00286 return false; 00287 } 00288 bool CbIsPointInGeom13(Vector3D const & P, int & ClrIdx) 00289 { 00290 if (P(0)>=0.0 && P(0)<=15.0 && P(1)>0.0 && P(1)<=10.0) return true; 00291 return false; 00292 } 00293 bool CbIsPointInGeom14(Vector3D const & P, int & ClrIdx) 00294 { 00295 // ground 00296 double L = 30.0; 00297 double H = 10.0; 00298 if (P(0)>=0.0 && P(0)<=L && P(1)>0.0 && P(1)<=H) 00299 { 00300 ClrIdx = -((static_cast<int>(P(1)) % 2) + 1); 00301 return true; 00302 } 00303 00304 // comet 00305 double xc = L/2.0; 00306 double yc = H+3.0; 00307 double r = 2.0; 00308 if (sqrt(pow(P(0)-xc,2.)+pow(P(1)-yc,2.))<=r) 00309 { 00310 ClrIdx = -3; 00311 return true; 00312 } 00313 00314 // default 00315 return false; 00316 } 00317 bool CbIsPointInGeom15(Vector3D const & P, int & ClrIdx) 00318 { 00319 double xc = 5.0; 00320 double yc = 5.0; 00321 double r = 4.0; 00322 if (sqrt(pow(P(0)-xc,2.)+pow(P(1)-yc,2.))<=r) return true; 00323 return false; 00324 } 00325 bool CbIsPointInGeom16(Vector3D const & P, int & ClrIdx) 00326 { 00327 if ((P(0)>0.5 && P(0)<1.5) && (P(1)>8.5 && P(1)<9.5)) return true; 00328 return false; 00329 } 00330 bool CbIsPointInGeom17(Vector3D const & P, int & ClrIdx) 00331 { 00332 if ((P(0)>0.5 && P(0)<1.5) && (P(1)>8.5 && P(1)<9.5)) 00333 { 00334 double xmin = 0.5; 00335 double l = 0.125; 00336 ClrIdx = -((static_cast<int>((P(0)-xmin)/l) % 2) + 1); 00337 return true; 00338 } 00339 return false; 00340 } 00341 00342 // Density callbacks 00343 double CbDensity1 (Vector3D const & P) { return 1.0; } 00344 double CbDensity2 (Vector3D const & P) { return 1.0; } 00345 double CbDensity3 (Vector3D const & P) { return 1.0; } 00346 double CbDensity4 (Vector3D const & P) { return 1.0; } 00347 double CbDensity5 (Vector3D const & P) { return 2.71; } // picogram/(micrometer*nanoseconds) 00348 double CbDensity6 (Vector3D const & P) { return 1.0; } 00349 double CbDensity7 (Vector3D const & P) { return 1.0; } 00350 double CbDensity8 (Vector3D const & P) { return 1.0; } 00351 double CbDensity10(Vector3D const & P) { return 1.0; } 00352 double CbDensity11(Vector3D const & P) { return 1.0; } 00353 double CbDensity12(Vector3D const & P) { return 1.0; } 00354 double CbDensity13(Vector3D const & P) { return 1.0; } 00355 double CbDensity14(Vector3D const & P) { return 1.0; } 00356 double CbDensity15(Vector3D const & P) { return 1.0; } 00357 double CbDensity16(Vector3D const & P) { return 1050.0; } // kg/m3 00358 double CbDensity17(Vector3D const & P) { return 1050.0; } // kg/m3 00359 00360 const double DTRUE = 1.0; 00361 00362 // Allocate model callbacks 00363 void CbAllocMdl1 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00364 { 00365 Prms.Set ("E nu psa", 4.0*PI*PI, 0.0, DTRUE); 00366 Inis.Set ("zero", DTRUE); 00367 Name = "NeoHookean"; 00368 } 00369 void CbAllocMdl2 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00370 { 00371 Prms.Set ("E nu psa", 100.0, 0.25, DTRUE); 00372 Inis.Set ("zero", DTRUE); 00373 Name = "NeoHookean"; 00374 } 00375 void CbAllocMdl3 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00376 { 00377 Prms.Set ("E nu psa", 1.0, 0.2, DTRUE); 00378 Inis.Set ("zero", DTRUE); 00379 Name = "NeoHookean"; 00380 } 00381 void CbAllocMdl4 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00382 { 00383 Prms.Set ("E nu psa", 210.0, 0.3, DTRUE); 00384 Inis.Set ("zero", DTRUE); 00385 Name = "NeoHookean"; 00386 } 00387 void CbAllocMdl5 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00388 { 00389 double E = 175.8; // E picogram/(micrometer*nanoseconds) == GPa 00390 double nu = 0.28; // nu 00391 Prms.Set ("E nu psa", E, nu, DTRUE); 00392 Inis.Set ("zero", DTRUE); 00393 Name = "NeoHookean"; 00394 } 00395 void CbAllocMdl6 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00396 { 00397 double xc=10.0; double yc=10.0; 00398 double d = sqrt(pow(P(0)-xc,2.0)+pow(P(1)-yc,2.0)); 00399 00400 if (d<=2.5) 00401 { 00402 Prms.Set ("E nu psa", 1.0, 0.49, DTRUE); 00403 Inis.Set ("zero", DTRUE); 00404 Name = "NeoHookean"; 00405 } 00406 else 00407 { 00408 Prms.Set ("E nu sY VM psa", 10000.0, 0.1, 2.0, DTRUE, DTRUE); 00409 Inis.Set ("zero", DTRUE); 00410 Name = "ElastoPlastic"; 00411 } 00412 } 00413 void CbAllocMdl7 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00414 { 00415 Prms.Set ("E nu psa", 1000.0, 0.25, DTRUE); 00416 Inis.Set ("zero", DTRUE); 00417 Name = "NeoHookean"; 00418 } 00419 void CbAllocMdl7b (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00420 { 00421 Prms.Set ("E nu sY VM psa", 100.0, 0.25, 3.0, DTRUE, DTRUE); 00422 Inis.Set ("zero", DTRUE); 00423 Name = "ElastoPlastic"; 00424 } 00425 void CbAllocMdl8 (Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00426 { 00427 Prms.Set ("E nu psa", 100.0, 0.25, DTRUE); 00428 Inis.Set ("zero", DTRUE); 00429 Name = "NeoHookean"; 00430 } 00431 void CbAllocMdl10(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00432 { 00433 Prms.Set ("E nu psa", 100.0, 0.0, DTRUE); 00434 Inis.Set ("zero", DTRUE); 00435 Name = "NeoHookean"; 00436 } 00437 void CbAllocMdl11(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00438 { 00439 Prms.Set ("E nu psa", 6778.0, 0.2103, DTRUE); 00440 Inis.Set ("sx sy sz", -30.0, -30.0, -30.0); 00441 Name = "NeoHookean"; 00442 } 00443 void CbAllocMdl12(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00444 { 00445 //Prms.Set ("E nu psa", 10.0, 0.3, DTRUE); 00446 //Inis.Set ("zero", DTRUE); 00447 //Name = "NeoHookean"; 00448 00449 Inis.Set ("zero", DTRUE); 00450 Prms.Set ("E nu sY VM psa", 100000.0, 0.3, 2.0, DTRUE, DTRUE); 00451 //Prms.Set ("E nu sY VM psa", 100.0, 0.3, 2.0, DTRUE, DTRUE); 00452 Name = "ElastoPlastic"; 00453 } 00454 void CbAllocMdl13(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00455 { 00456 Prms.Set ("E nu psa", 10.0, 0.3, DTRUE); 00457 Inis.Set ("zero", DTRUE); 00458 Name = "NeoHookean"; 00459 } 00460 void CbAllocMdl14(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00461 { 00462 // comet 00463 double L = 30.0; 00464 double H = 10.0; 00465 double xc = L/2.0; 00466 double yc = H+3.0; 00467 double r = 2.0; 00468 if (sqrt(pow(P(0)-xc,2.)+pow(P(1)-yc,2.))<=r) 00469 { 00470 Prms.Set ("E nu psa", 100000.0, 0.1, DTRUE); 00471 Inis.Set ("zero", DTRUE); 00472 Name = "NeoHookean"; 00473 Tag = -1; 00474 } 00475 00476 // ground 00477 else 00478 { 00479 Inis.Set ("zero", DTRUE); 00480 Prms.Set ("E nu sY VM psa", 100.0, 0.3, 2.0, DTRUE, DTRUE); 00481 //Prms.Set ("E nu c phi DP psa", 100.0, 0.3, 0.0, 20.0, DTRUE, DTRUE); 00482 Name = "ElastoPlastic"; 00483 Tag = -2; 00484 } 00485 } 00486 void CbAllocMdl15(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00487 { 00488 Prms.Set ("E nu psa", 1000.0, 0.1, DTRUE); 00489 Inis.Set ("zero", DTRUE); 00490 Name = "NeoHookean"; 00491 } 00492 void CbAllocMdl16(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00493 { 00494 Prms.Set ("E nu", 1.0e+6, 0.3, DTRUE); 00495 Inis.Set ("zero", DTRUE); 00496 Name = "NeoHookean"; 00497 } 00498 void CbAllocMdl17(Vector3D const & P, String & Name, int & Tag, int NDim, SDPair & Prms, SDPair & Inis, String & Scheme) 00499 { 00500 Prms.Set ("E nu", 1.0e+6, 0.3, DTRUE); 00501 Inis.Set ("zero", DTRUE); 00502 Name = "NeoHookean"; 00503 } 00504 00505 // Initial velocity 00506 void CbIniVeloc1 (Vector3D const & P, Vector3D & v) { v = 0.1, 0.0, 0.0; } 00507 void CbIniVeloc2 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00508 void CbIniVeloc3 (Vector3D const & P, Vector3D & v) 00509 { 00510 /* 00511 if (P(0)<0.5 && P(1)<0.5 && P(0)>0.1 && P(1)>0.1) v = 0.15, 0.15, 0.0; 00512 else if (P(0)>0.5 && P(1)>0.5 && P(0)<0.9 && P(1)<0.9) v = -0.15, -0.1, 0.0; 00513 else v = 0.0, 0.0, 0.0; 00514 */ 00515 /* 00516 if (P(0)<0.5 && P(1)<0.5 && P(0)>0.1 && P(1)>0.1) v = 0.1, 0.1, 0.0; 00517 else if (P(0)>0.5 && P(1)>0.5 && P(0)<0.9 && P(1)<0.9) v = -0.1, -0.1, 0.0; 00518 else v = 0.0, 0.0, 0.0; 00519 */ 00520 if (P(0)<0.325 && P(1)<0.325) v = 0.1, 0.1, 0.0; 00521 else v = -0.1, -0.1, 0.0; 00522 } 00523 void CbIniVeloc4 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00524 void CbIniVeloc5 (Vector3D const & P, Vector3D & v) { v = 0.0; } // micrometer/nanoseconds 00525 void CbIniVeloc6 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00526 void CbIniVeloc7 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00527 void CbIniVeloc8 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00528 void CbIniVeloc10(Vector3D const & P, Vector3D & v) 00529 { 00530 double v0 = 0.1; // initial velocity 00531 //double E = 100; // Young modulus 00532 //double rho0 = 1.0; // initial density 00533 double mode = 5.0; // Mode 00534 double bet = ((2.0*mode-1.0)/2.0)*PI/25.0; // Eigenvalues 00535 //double c = sqrt(E/rho0); // elastic wave speed 00536 //double w = bet*c; // frequency 00537 double vel0 = v0*sin(bet*P(0)); // initial velocity 00538 v = vel0, 0.0, 0.0; 00539 //double xcm = 12.5; // centre of mass 00540 //double velcm= function(t) v0*cos(w*t)*sin(bet*xcm) // correct velocity (center of mass) 00541 } 00542 void CbIniVeloc11 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00543 void CbIniVeloc12 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00544 void CbIniVeloc13 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00545 void CbIniVeloc14 (Vector3D const & P, Vector3D & v) 00546 { 00547 // comet 00548 double L = 30.0; 00549 double H = 10.0; 00550 double xc = L/2.0; 00551 double yc = H+3.0; 00552 double r = 2.0; 00553 if (sqrt(pow(P(0)-xc,2.)+pow(P(1)-yc,2.))<=r) v = 0.0, -10.0, 0.0; 00554 else v = 0.0; 00555 } 00556 void CbIniVeloc15 (Vector3D const & P, Vector3D & v) 00557 { 00558 // comet 00559 double xc = 5.0; 00560 double yc = 5.0; 00561 double dfdx = 2.0*(P(0)-xc); 00562 double dfdy = 2.0*(P(1)-yc); 00563 double d = sqrt(pow(P(0)-xc,2.)+pow(P(1)-yc,2.)); 00564 v = -dfdy, dfdx, 0.0; 00565 v *= (0.1*d*NormV(v)); 00566 } 00567 void CbIniVeloc16 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00568 void CbIniVeloc17 (Vector3D const & P, Vector3D & v) { v = 0.0; } 00569 00570 /* x sqrt(x) 2*sqrt(x) log2(2*sqrt(x)) 1/sqrt(x) 1-1/sqrt(x) 00571 * 1 1 2 1 1 0.0 00572 * 4 2 4 2 0.5 0.5 00573 * 16 4 8 3 0.25 0.75 00574 * 00575 * 1 4 16 00576 * 1_ _______ 1_ _______ 1_ _______ 00577 * | | 0.5_| | _ |_3L/4 0.75>|_|_|_|_|<7L/8 00578 * 0_| _ |_L/2 0_|___|___| 0_|_|_|_|_| 00579 * | | | | _ |_L/4 |_|_|_|_| 00580 * -1_|_______| -1_|___|___| -1_|_|_|_|_|<L/8 */ 00581 00582 // Applied tractions 00583 bool CbHasTraction1 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00584 bool CbHasTraction2 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00585 { 00586 double xmax = 1.0; 00587 double ymax = 1.0; 00588 if (N[0](1)>ymax-1.01*L(1) && N[1](0)<xmax/1.99) // half/upper cells 00589 { 00590 double sqnpc = sqrt(static_cast<double>(nPCell)); 00591 if (P(1)>N[3](1)-0.51*L(1)/sqnpc) 00592 { 00593 double q = 5.0; 00594 t = 0.0, -q*L(0)/sqnpc, 0.0; 00595 return true; 00596 } 00597 } 00598 return false; 00599 } 00600 bool CbHasTraction3 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00601 bool CbHasTraction4 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00602 { 00603 if (N[0](1)>0.99-L(1)) 00604 { 00605 if (P(1)>N[0](1)+L(1)/2.01) 00606 { 00607 double p = 1.0; 00608 if (nPCell==1) t = 0.0, p*L(0) , 0.0; 00609 else t = 0.0, p*L(0)/2.0, 0.0; 00610 return true; 00611 } 00612 } 00613 return false; 00614 } 00615 bool CbHasTraction5 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00616 { 00617 if (N[0](0)>59.9-L(0)) 00618 { 00619 if (P(0)>59.9-L(0)/2.0) 00620 { 00621 double p = 400.0; // 40 picogram/(micrometer*nanoseconds) == GPa 00622 double val = p*L(1)/sqrt(nPCell); 00623 t = val, 0.0, 0.0; 00624 return true; 00625 } 00626 } 00627 return false; 00628 } 00629 bool CbHasTraction6 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00630 { 00631 double xc=10.0; double yc=10.0; 00632 double f = 1000.0; 00633 double d = sqrt(pow(P(0)-xc,2.0)+pow(P(1)-yc,2.0)); 00634 00635 if (d<=2.5) 00636 { 00637 Vector3D c, r; 00638 c = xc, yc, 0.0; 00639 r = P - c; 00640 double n = sqrt(dot(r,r)); 00641 if (n>0.0) r /= n; 00642 else r = 0.0; 00643 t = f*r; 00644 return true; 00645 } 00646 else return false; 00647 } 00648 bool CbHasTraction7 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00649 { 00650 double npc = sqrt(static_cast<double>(nPCell)); 00651 if (P(0)>19.0-L(0)/npc) 00652 { 00653 t = 0.0, -10.0, 0.0; 00654 return true; 00655 } 00656 return false; 00657 } 00658 bool CbHasTraction8 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00659 { 00660 double npc = sqrt(static_cast<double>(nPCell)); 00661 if (P(0)>5.0 && P(0)<5.0+L(0)/npc && P(1)>12.0-L(1)/npc) 00662 { 00663 t = 0.0, -10.0, 0.0; 00664 return true; 00665 } 00666 return false; 00667 } 00668 bool CbHasTraction10(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00669 bool CbHasTraction11(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00670 { 00671 double xmax = 10.0; 00672 double ymax = 10.0; 00673 double sqnpc = sqrt(static_cast<double>(nPCell)); 00674 if (N[0](1)>ymax-1.01*L(1) && 0.99*N[1](0)<xmax) // upper cells 00675 { 00676 if (P(1)>N[3](1)-0.51*L(1)/sqnpc) 00677 { 00678 double q = 30.0; 00679 t = 0.0, -q*L(0)/sqnpc, 0.0; 00680 return true; 00681 } 00682 } 00683 if (N[0](0)>xmax-1.01*L(0) && 0.99*N[1](1)<ymax) // right cells 00684 { 00685 if (P(0)>N[1](0)-0.51*L(0)/sqnpc) 00686 { 00687 double q = 30.0; 00688 t = -q*L(0)/sqnpc, 0.0, 0.0; 00689 return true; 00690 } 00691 } 00692 return false; 00693 } 00694 bool CbHasTraction12(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00695 bool CbHasTraction13(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00696 bool CbHasTraction14(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00697 bool CbHasTraction15(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00698 bool CbHasTraction16(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00699 bool CbHasTraction17(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00700 00701 // Applied displacements 00702 bool CbHasAppDisp1 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00703 bool CbHasAppDisp2 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00704 bool CbHasAppDisp3 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00705 bool CbHasAppDisp4 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00706 bool CbHasAppDisp5 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00707 bool CbHasAppDisp6 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00708 bool CbHasAppDisp7 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00709 bool CbHasAppDisp8 (Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00710 bool CbHasAppDisp10(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00711 bool CbHasAppDisp11(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00712 bool CbHasAppDisp12(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00713 bool CbHasAppDisp13(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) 00714 { 00715 if (P(0)>=0.0 && P(0)<=0.5 && P(1)>9.7) 00716 { 00717 t = 0., -0.1, 0.; // downward velocity 00718 return true; 00719 } 00720 return false; 00721 } 00722 bool CbHasAppDisp14(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00723 bool CbHasAppDisp15(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00724 bool CbHasAppDisp16(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00725 bool CbHasAppDisp17(Vector3D const & P, Array<Vector3D> const & N, Vector3D const & L, int nPCell, Vector3D & t) { return false; } 00726 00727 // Body mass callbacks 00728 static inline void CbB1 (double t, Vector3D & B) { B=0.0; } 00729 static inline void CbB2 (double t, Vector3D & B) { B=0.0; } 00730 static inline void CbB3 (double t, Vector3D & B) { B=0.0; } 00731 static inline void CbB4 (double t, Vector3D & B) { B=0.0; } 00732 static inline void CbB5 (double t, Vector3D & B) { B=0.0; } 00733 static inline void CbB6 (double t, Vector3D & B) { B=0.0; }//, -1.0, 0.0; } 00734 static inline void CbB7 (double t, Vector3D & B) { B=0.0; }//, -1.0, 0.0; } 00735 static inline void CbB8 (double t, Vector3D & B) { B=0.0; }//, -1.0, 0.0; } 00736 static inline void CbB10(double t, Vector3D & B) { B=0.0; } 00737 static inline void CbB11(double t, Vector3D & B) { B=0.0; } 00738 static inline void CbB12(double t, Vector3D & B) 00739 { 00740 //double tsw = 5.0; // very funny! 00741 double tsw = 20.0; 00742 double mul = (t<tsw ? sin(PI*t/(2.0*tsw)) : 1.0); 00743 B = 0.0, -mul, 0.0; 00744 } 00745 static inline void CbB13(double t, Vector3D & B) { B=0.0; } 00746 static inline void CbB14(double t, Vector3D & B) { B=0.0; } 00747 static inline void CbB15(double t, Vector3D & B) { B=0.0; } 00748 static inline void CbB16(double t, Vector3D & B) { B=0.0, -1000.0, 0.0; } // m/s2 00749 static inline void CbB17(double t, Vector3D & B) { B=0.0, -3500.0, 0.0; } // m/s2 00750 00751 // Loading multiplier callbacks 00752 static inline void CbLdM1 (double t, double & M) { M=0.0; } // t in seconds (s) 00753 static inline void CbLdM2 (double t, double & M) { M=t; } // t in seconds (s) 00754 static inline void CbLdM3 (double t, double & M) { M=0.0; } // t in seconds (s) 00755 static inline void CbLdM4 (double t, double & M) { M=(t<0.5 ? t/0.5 : 1.0); } // t in seconds (s) 00756 static inline void CbLdM5 (double t, double & M) { M=(t<60.0 ? t/60.0 : 1.0); } // t in nano-seconds (ns) 00757 //static inline void CbLdM5 (double t, double & M) { M=(t<10.0 ? t/10.0 : 1.0); } // t in nano-seconds (ns) 00758 //static inline void CbLdM6 (double t, double & M) { M=(t>0.5 ? (t<1 ? 1.0 : 0.0) : 0.0); } // t in seconds (s) 00759 static inline void CbLdM6 (double t, double & M) { M=1.0; }//(t>0.01 ? (t<0.11 ? 1.0 : 0.0) : 0.0); } // t in seconds (s) 00760 static inline void CbLdM7 (double t, double & M) { M=(t<10.0 ? t/10.0 : 1.0); } // t in seconds (s) 00761 static inline void CbLdM8 (double t, double & M) { M=(t<5.0 ? t/5.0 : 0.0); } // t in seconds (s) 00762 static inline void CbLdM10(double t, double & M) { M=0.0; } // t in seconds (s) 00763 static inline void CbLdM11(double t, double & M) { M=(t<100.0 ? t/100.0 : 1.0); } // t in seconds (s) 00764 static inline void CbLdM12(double t, double & M) { M=1.0; } // t in seconds (s) 00765 static inline void CbLdM13(double t, double & M) { M=(t<10.0 ? t/10.0 : 1.0); } // t in seconds (s) 00766 static inline void CbLdM14(double t, double & M) { M=0.0; } // t in seconds (s) 00767 static inline void CbLdM15(double t, double & M) { M=0.0; } // t in seconds (s) 00768 static inline void CbLdM16(double t, double & M) { M=0.0; } // t in seconds (s) 00769 static inline void CbLdM17(double t, double & M) { M=0.0; } // t in seconds (s) 00770 00771 // Correct velocities callbacks 00772 static inline void CbVel1 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.1*cos(2.0*PI*t),0,0; } 00773 static inline void CbVel2 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00774 static inline void CbVel3 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00775 static inline void CbVel4 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00776 static inline void CbVel5 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00777 static inline void CbVel6 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00778 static inline void CbVel7 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00779 static inline void CbVel8 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00780 static inline void CbVel10(double t, Vector3D const & XY, Vector3D & Vel) 00781 { 00782 double v0 = 0.1; // initial velocity 00783 double E = 100; // Young modulus 00784 double rho0 = 1.0; // initial density 00785 double mode = 5.0; // Mode 00786 double bet = ((2.0*mode-1.0)/2.0)*PI/25.0; // Eigenvalues 00787 double c = sqrt(E/rho0); // elastic wave speed 00788 double w = bet*c; // frequency 00789 double xcm = 12.5; // centre of mass 00790 double velcm = v0*cos(w*t)*sin(bet*xcm); // correct velocity (center of mass) 00791 Vel = velcm, 0.0, 0.0; 00792 } 00793 static inline void CbVel11 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00794 static inline void CbVel12 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00795 static inline void CbVel13 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00796 static inline void CbVel14 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00797 static inline void CbVel15 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00798 static inline void CbVel16 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00799 static inline void CbVel17 (double t, Vector3D const & XY, Vector3D & Vel) { Vel=0.0; } 00800 00801 // Correct stresses callbacks 00802 static inline void CbStress1 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00803 static inline void CbStress2 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00804 static inline void CbStress3 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00805 static inline void CbStress4 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00806 static inline void CbStress5 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00807 static inline void CbStress6 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00808 static inline void CbStress7 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00809 static inline void CbStress8 (double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00810 static inline void CbStress10(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00811 static inline void CbStress11(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00812 static inline void CbStress12(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00813 static inline void CbStress13(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00814 static inline void CbStress14(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00815 static inline void CbStress15(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00816 static inline void CbStress16(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00817 static inline void CbStress17(double t, Vector3D const & XY, STensor2 & Sig) { Sig=0.0; } 00818 00820 struct Problem2D 00821 { 00822 // Data 00823 int Id; 00824 int NPcell; 00825 Grid2D * G; 00826 MPoints2D * P; 00827 ptIsPointInGeom pGeom; 00828 ptB pB; 00829 ptLdM pM; 00830 ptVeloc pVeloc; 00831 ptStress pStress; 00832 bool CPDI; 00833 bool MPM; 00834 bool USF; 00835 bool SmallSt; 00836 bool FwEuler; 00837 bool SaveVTU; 00838 bool SavePNG; 00839 bool SavePointData; 00840 double t; 00841 double tf; 00842 double Dt; 00843 double Dtout; 00844 size_t Outp; 00845 Vector3D B; 00846 double M; 00847 00848 // Output arrays 00849 Array<double> Out_t; // time array 00850 Array<double> Out_Wx_t; // Body weights (in time) 00851 Array<double> Out_sE_t; // Out strain Energy (in time) 00852 Array<double> Out_kE_t; // Out kinetic Energy (in time) 00853 Array<double> Out_tE_t; // Out total Energy (in time) 00854 Array<Array<double> > Out_p_vx_t; // Out particle x-velocity (in time) Size=np 00855 Array<Array<double> > Out_p_vy_t; // Out particle y-velocity (in time) Size=np 00856 Array<Array<double> > Out_p_cvx_t; // Out correct particle x-velocity (in time) Size=np 00857 Array<Array<double> > Out_p_sxx_t; // Out particle xx-stress (in time) Size=np 00858 Array<Array<double> > Out_p_syy_t; // Out particle yy-stress (in time) Size=np 00859 Array<Array<double> > Out_p_exx_t; // Out particle xx-strain (in time) Size=np 00860 Array<Array<double> > Out_p_eyy_t; // Out particle yy-strain (in time) Size=np 00861 Array<Array<double> > Out_p_fx_t; // Out external force x (in time) Size=np 00862 Array<Array<double> > Out_p_fy_t; // Out external force y (in time) Size=np 00863 Array<Array<double> > Out_p_ux_t; // Out displacement x (in time) Size=np 00864 Array<Array<double> > Out_p_uy_t; // Out displacement y (in time) Size=np 00865 }; 00866 00868 void ReSetOutputArrays (Problem2D & Prob) 00869 { 00870 // Erase output arrays 00871 Prob.Out_t .Resize(0); 00872 Prob.Out_Wx_t .Resize(0); 00873 Prob.Out_sE_t .Resize(0); 00874 Prob.Out_kE_t .Resize(0); 00875 Prob.Out_tE_t .Resize(0); 00876 Prob.Out_p_vx_t .Resize(Prob.P->nPoints()); 00877 Prob.Out_p_vy_t .Resize(Prob.P->nPoints()); 00878 Prob.Out_p_cvx_t.Resize(Prob.P->nPoints()); 00879 Prob.Out_p_sxx_t.Resize(Prob.P->nPoints()); 00880 Prob.Out_p_syy_t.Resize(Prob.P->nPoints()); 00881 Prob.Out_p_exx_t.Resize(Prob.P->nPoints()); 00882 Prob.Out_p_eyy_t.Resize(Prob.P->nPoints()); 00883 Prob.Out_p_fx_t .Resize(Prob.P->nPoints()); 00884 Prob.Out_p_fy_t .Resize(Prob.P->nPoints()); 00885 Prob.Out_p_ux_t .Resize(Prob.P->nPoints()); 00886 Prob.Out_p_uy_t .Resize(Prob.P->nPoints()); 00887 } 00888 00889 Fl_Menu_Item Problems[] = 00890 { 00891 {"1. Linear vibration" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00892 {"2. Indentation (WCCM8)" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00893 {"3. Bouncing disks" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00894 {"4. 2D vibration - small" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00895 {"5. Particle separation problem" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00896 {"6. Explosion" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00897 {"7. Cantilever" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00898 {"8. Column (buckling)" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00899 {"9. Cantilever (plastic)" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00900 {"10.Linear vibration (continuum)" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00901 {"11.Cylindrical Hole" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00902 {"12.Spoil pile" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00903 {"13.Footing penetration" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00904 {"14.Comet impact" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00905 {"15.Rotating disk" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00906 {"16. SBB-02a: gravity" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00907 {"17. SBB-02b: gravity" , 0, 0, 0, 0, FL_NORMAL_LABEL, 0, 0, 0}, 00908 {0,0,0,0,0,0,0,0,0} 00909 }; 00910 size_t DefaultProblem = 3; 00911 00913 void SetProblem (Problem2D & Prob) 00914 { 00915 // Deallocate previous grid and matpoints 00916 if (Prob.G!=NULL) delete Prob.G; 00917 if (Prob.P!=NULL) delete Prob.P; 00918 00919 // Set problem data 00920 switch (Prob.Id) 00921 { 00922 case 1: 00923 { 00924 // Allocate grid 00925 Prob.G = new Grid2D (/*xMin*/0, /*yMin*/0, /*nRow*/4, /*nCol*/4, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed1); 00926 00927 // Allocate material points 00928 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom1, &CbDensity1, &CbAllocMdl1, &CbIniVeloc1, &CbHasTraction1, &CbHasAppDisp1); 00929 00930 // Set callbacks 00931 Prob.pGeom = &CbIsPointInGeom1; 00932 Prob.pB = &CbB1; 00933 Prob.pM = &CbLdM1; 00934 Prob.pVeloc = &CbVel1; 00935 Prob.pStress = &CbStress1; 00936 00937 // Set constants 00938 Prob.t = 0.0; // current time 00939 Prob.tf = 5.0; // final time 00940 Prob.Dt = 0.015;// 0.1/(2.0*PI); // time increment 00941 Prob.Dtout = 0.015; // time increment for output 00942 Prob.Outp = 0; // point/particle number to output data 00943 Prob.B = 0.0; // current b-body mass 00944 Prob.M = 0.0; // multiplier for external forces 00945 00946 break; 00947 } 00948 case 2: 00949 { 00950 // Allocate grid 00951 Prob.G = new Grid2D (/*xMin*/-0.25, /*yMin*/-0.25, /*nRow*/7, /*nCol*/7, /*Lx*/0.25, /*Ly*/0.25, &CbIsFixed2); 00952 00953 // Allocate material points 00954 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom2, &CbDensity2, &CbAllocMdl2, &CbIniVeloc2, &CbHasTraction2, &CbHasAppDisp2); 00955 00956 // Set callbacks 00957 Prob.pGeom = &CbIsPointInGeom2; 00958 Prob.pB = &CbB2; 00959 Prob.pM = &CbLdM2; 00960 Prob.pVeloc = &CbVel2; 00961 Prob.pStress = &CbStress2; 00962 00963 // Set constants 00964 Prob.t = 0.0; // current time 00965 Prob.tf = 8.0; // final time 00966 Prob.Dt = 0.01; // time increment 00967 Prob.Dtout = 0.1; // time increment for output 00968 Prob.Outp = 0; // point/particle number to output data 00969 Prob.B = 0.0; // current b-body mass 00970 Prob.M = 0.0; // multiplier for external forces 00971 00972 break; 00973 } 00974 case 3: 00975 { 00976 // Allocate grid 00977 //Prob.G = new Grid2D (/*xMin*/-0.10, /*yMin*/-0.10, /*nRow*/24, /*nCol*/24, /*Lx*/0.05, /*Ly*/0.05, &CbIsFixed3); 00978 Prob.G = new Grid2D (/*xMin*/-0.10, /*yMin*/-0.10, /*nRow*/18, /*nCol*/18, /*Lx*/0.05, /*Ly*/0.05, &CbIsFixed3); 00979 00980 // Allocate material points 00981 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom3, &CbDensity3, &CbAllocMdl3, &CbIniVeloc3, &CbHasTraction3, &CbHasAppDisp3); 00982 00983 // Set callbacks 00984 Prob.pGeom = &CbIsPointInGeom3; 00985 Prob.pB = &CbB3; 00986 Prob.pM = &CbLdM3; 00987 Prob.pVeloc = &CbVel3; 00988 Prob.pStress = &CbStress3; 00989 00990 // Set constants 00991 Prob.t = 0.0; // current time 00992 Prob.tf = 10.0; // final time 00993 Prob.Dt = 0.001; // time increment 00994 Prob.Dtout = 0.1; // time increment for output 00995 Prob.Outp = 0; // point/particle number to output data 00996 Prob.B = 0.0; // current b-body mass 00997 Prob.M = 0.0; // multiplier for external forces 00998 00999 /* 01000 // Boundary conditions 01001 for (int i=0; i<Prob.G->nNodes(); ++i) 01002 { 01003 if (Prob.G->N(i)(0)<0.10) Prob.G->SetFixed(i,FIX_XY); // left 01004 if (Prob.G->N(i)(0)>0.89) Prob.G->SetFixed(i,FIX_XY); // right 01005 if (Prob.G->N(i)(1)<0.10) Prob.G->SetFixed(i,FIX_XY); // bottom 01006 if (Prob.G->N(i)(1)>0.89) Prob.G->SetFixed(i,FIX_XY); // top 01007 } 01008 */ 01009 01010 break; 01011 } 01012 case 4: 01013 { 01014 // Allocate grid 01015 Prob.G = new Grid2D (/*xMin*/-0.25, /*yMin*/-0.25, /*nRow*/7, /*nCol*/7, /*Lx*/0.25, /*Ly*/0.25, &CbIsFixed4); 01016 01017 // Allocate material points 01018 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom4, &CbDensity4, &CbAllocMdl4, &CbIniVeloc4, &CbHasTraction4, &CbHasAppDisp4); 01019 01020 // Set callbacks 01021 Prob.pGeom = &CbIsPointInGeom4; 01022 Prob.pB = &CbB4; 01023 Prob.pM = &CbLdM4; 01024 Prob.pVeloc = &CbVel4; 01025 Prob.pStress = &CbStress4; 01026 01027 // Set constants 01028 Prob.t = 0.0; // current time 01029 Prob.tf = 4.0; // final time 01030 Prob.Dt = 0.0005; // time increment 01031 Prob.Dtout = 0.01; // time increment for output 01032 Prob.Outp = 0; // point/particle number to output data 01033 Prob.B = 0.0; // current b-body mass 01034 Prob.M = 0.0; // multiplier for external forces 01035 01036 break; 01037 } 01038 case 5: 01039 { 01040 // Allocate grid 01041 Prob.G = new Grid2D (/*xMin*/-2.0, /*yMin*/-2.0, /*nRow*/23, /*nCol*/123, /*Lx*/2.0, /*Ly*/2.0, &CbIsFixed5); 01042 01043 // Allocate material points 01044 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom5, &CbDensity5, &CbAllocMdl5, &CbIniVeloc5, &CbHasTraction5, &CbHasAppDisp5); 01045 01046 // Set callbacks 01047 Prob.pGeom = &CbIsPointInGeom5; 01048 Prob.pB = &CbB5; 01049 Prob.pM = &CbLdM5; 01050 Prob.pVeloc = &CbVel5; 01051 Prob.pStress = &CbStress5; 01052 01053 // Set constants 01054 Prob.t = 0.0; // current time (nano-seconds) 01055 Prob.tf = 60.0; // final time (nano-seconds) 01056 Prob.Dt = 0.02; // time increment (nano-seconds) 01057 Prob.Dtout = 0.1; // time increment for output (nano-seconds) 01058 Prob.Outp = 0; // point/particle number to output data 01059 Prob.B = 0.0; // current b-body mass 01060 Prob.M = 0.0; // multiplier for external forces 01061 01062 break; 01063 } 01064 case 6: 01065 { 01066 // Allocate grid 01067 Prob.G = new Grid2D (/*xMin*/-1.0, /*yMin*/-1.0, /*nRow*/23, /*nCol*/23, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed6); 01068 01069 // Allocate material points 01070 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom6, &CbDensity6, &CbAllocMdl6, &CbIniVeloc6, &CbHasTraction6, &CbHasAppDisp6); 01071 01072 // Set callbacks 01073 Prob.pGeom = &CbIsPointInGeom6; 01074 Prob.pB = &CbB6; 01075 Prob.pM = &CbLdM6; 01076 Prob.pVeloc = &CbVel6; 01077 Prob.pStress = &CbStress6; 01078 01079 // Set constants 01080 Prob.t = 0.0; // current time 01081 Prob.tf = 1.0; // final time 01082 Prob.Dt = 0.0001; // time increment 01083 Prob.Dtout = 0.01; // time increment for output 01084 Prob.Outp = 0; // point/particle number to output data 01085 Prob.B = 0.0; // current b-body mass 01086 Prob.M = 0.0; // multiplier for external forces 01087 01088 break; 01089 } 01090 case 7: 01091 { 01092 // Allocate grid 01093 Prob.G = new Grid2D (/*xMin*/-1.0, /*yMin*/-1.0, /*nRow*/16, /*nCol*/23, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed7); 01094 01095 // Allocate material points 01096 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom7, &CbDensity7, &CbAllocMdl7, &CbIniVeloc7, &CbHasTraction7, &CbHasAppDisp7); 01097 01098 // Set callbacks 01099 Prob.pGeom = &CbIsPointInGeom7; 01100 Prob.pB = &CbB7; 01101 Prob.pM = &CbLdM7; 01102 Prob.pVeloc = &CbVel7; 01103 Prob.pStress = &CbStress7; 01104 01105 // Set constants 01106 Prob.t = 0.0; // current time 01107 Prob.tf = 10.0; // final time 01108 Prob.Dt = 0.001; // time increment 01109 Prob.Dtout = 0.1; // time increment for output 01110 Prob.Outp = 0; // point/particle number to output data 01111 Prob.B = 0.0; // current b-body mass 01112 Prob.M = 0.0; // multiplier for external forces 01113 01114 break; 01115 } 01116 case 8: 01117 { 01118 // Allocate grid 01119 Prob.G = new Grid2D (/*xMin*/-3.0, /*yMin*/-1.0, /*nRow*/16, /*nCol*/20, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed8); 01120 01121 // Allocate material points 01122 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom8, &CbDensity8, &CbAllocMdl8, &CbIniVeloc8, &CbHasTraction8, &CbHasAppDisp8); 01123 01124 // Set callbacks 01125 Prob.pGeom = &CbIsPointInGeom8; 01126 Prob.pB = &CbB8; 01127 Prob.pM = &CbLdM8; 01128 Prob.pVeloc = &CbVel8; 01129 Prob.pStress = &CbStress8; 01130 01131 // Set constants 01132 Prob.t = 0.0; // current time 01133 Prob.tf = 20.0; // final time 01134 Prob.Dt = 0.001; // time increment 01135 Prob.Dtout = 0.1; // time increment for output 01136 Prob.Outp = 0; // point/particle number to output data 01137 Prob.B = 0.0; // current b-body mass 01138 Prob.M = 0.0; // multiplier for external forces 01139 01140 break; 01141 } 01142 case 9: 01143 { 01144 // Allocate grid 01145 Prob.G = new Grid2D (/*xMin*/-1.0, /*yMin*/-1.0, /*nRow*/16, /*nCol*/23, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed7); 01146 01147 // Allocate material points 01148 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom7, &CbDensity7, &CbAllocMdl7b, &CbIniVeloc7, &CbHasTraction7, &CbHasAppDisp7); 01149 01150 // Set callbacks 01151 Prob.pGeom = &CbIsPointInGeom7; 01152 Prob.pB = &CbB7; 01153 Prob.pM = &CbLdM7; 01154 Prob.pVeloc = &CbVel7; 01155 Prob.pStress = &CbStress7; 01156 01157 // Set constants 01158 Prob.t = 0.0; // current time 01159 Prob.tf = 10.0; // final time 01160 Prob.Dt = 0.001; // time increment 01161 Prob.Dtout = 0.1; // time increment for output 01162 Prob.Outp = 0; // point/particle number to output data 01163 Prob.B = 0.0; // current b-body mass 01164 Prob.M = 0.0; // multiplier for external forces 01165 01166 break; 01167 } 01168 case 10: 01169 { 01170 // Allocate grid 01171 Prob.G = new Grid2D (/*xMin*/-1.0, /*yMin*/0, /*nRow*/4, /*nCol*/28, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed10); 01172 01173 // Allocate material points 01174 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom10, &CbDensity10, &CbAllocMdl10, &CbIniVeloc10, &CbHasTraction10, &CbHasAppDisp10); 01175 01176 // Set callbacks 01177 Prob.pGeom = &CbIsPointInGeom10; 01178 Prob.pB = &CbB10; 01179 Prob.pM = &CbLdM10; 01180 Prob.pVeloc = &CbVel10; 01181 Prob.pStress = &CbStress10; 01182 01183 // Set constants 01184 Prob.t = 0.0; // current time 01185 Prob.tf = 100.0; // final time 01186 Prob.Dt = 0.01; // time increment 01187 Prob.Dtout = 1.0; // time increment for output 01188 Prob.Outp = 12; // point/particle number to output data 01189 Prob.B = 0.0; // current b-body mass 01190 Prob.M = 0.0; // multiplier for external forces 01191 01192 break; 01193 } 01194 case 11: 01195 { 01196 /* Midas: GeoStrucAnal_01.pdf */ 01197 01198 // Allocate grid 01199 Prob.G = new Grid2D (/*xMin*/-1.0, /*yMin*/-1.0, /*nRow*/13, /*nCol*/13, /*Lx*/1.0, /*Ly*/1.0, &CbIsFixed11); 01200 01201 // Allocate material points 01202 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom11, &CbDensity11, &CbAllocMdl11, &CbIniVeloc11, &CbHasTraction11, &CbHasAppDisp11); 01203 01204 // Set callbacks 01205 Prob.pGeom = &CbIsPointInGeom11; 01206 Prob.pB = &CbB11; 01207 Prob.pM = &CbLdM11; 01208 Prob.pVeloc = &CbVel11; 01209 Prob.pStress = &CbStress11; 01210 01211 // Set constants 01212 Prob.t = 0.0; // current time 01213 Prob.tf = 100.0; // final time 01214 Prob.Dt = 0.01; // time increment 01215 Prob.Dtout = 1.0; // time increment for output 01216 Prob.Outp = 0; // point/particle number to output data 01217 Prob.B = 0.0; // current b-body mass 01218 Prob.M = 0.0; // multiplier for external forces 01219 01220 break; 01221 } 01222 case 12: 01223 { 01224 double H = 10.0; // height 01225 double L = 40.0; // length 01226 int ndivx = 40; 01227 int ndivy = 10; 01228 double Lx = L/ndivx; 01229 double Ly = H/ndivy; 01230 01231 // Allocate grid 01232 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed12); 01233 01234 // Allocate material points 01235 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom12, &CbDensity12, &CbAllocMdl12, &CbIniVeloc12, &CbHasTraction12, &CbHasAppDisp12); 01236 01237 // Set callbacks 01238 Prob.pGeom = &CbIsPointInGeom12; 01239 Prob.pB = &CbB12; 01240 Prob.pM = &CbLdM12; 01241 Prob.pVeloc = &CbVel12; 01242 Prob.pStress = &CbStress12; 01243 01244 // Set constants 01245 Prob.t = 0.0; // current time 01246 Prob.tf = 30.0; // final time 01247 Prob.Dt = 0.01; // time increment 01248 Prob.Dtout = 0.1; // time increment for output 01249 Prob.Outp = 0; // point/particle number to output data 01250 Prob.B = 0.0; // current b-body mass 01251 Prob.M = 0.0; // multiplier for external forces 01252 01253 break; 01254 } 01255 case 13: 01256 { 01257 double L = 15.0; // length 01258 double H = 10.0; // height 01259 int ndivx = 60; 01260 int ndivy = 40; 01261 double Lx = L/ndivx; 01262 double Ly = H/ndivy; 01263 01264 // Allocate grid 01265 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3+5, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed13); 01266 01267 // Allocate material points 01268 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom13, &CbDensity13, &CbAllocMdl13, &CbIniVeloc13, &CbHasTraction13, &CbHasAppDisp13); 01269 01270 // Set callbacks 01271 Prob.pGeom = &CbIsPointInGeom13; 01272 Prob.pB = &CbB13; 01273 Prob.pM = &CbLdM13; 01274 Prob.pVeloc = &CbVel13; 01275 Prob.pStress = &CbStress13; 01276 01277 // Set constants 01278 Prob.t = 0.0; // current time 01279 Prob.tf = 30.0; // final time 01280 Prob.Dt = 0.01; // time increment 01281 Prob.Dtout = 0.1; // time increment for output 01282 Prob.Outp = 0; // point/particle number to output data 01283 Prob.B = 0.0; // current b-body mass 01284 Prob.M = 0.0; // multiplier for external forces 01285 01286 break; 01287 } 01288 case 14: 01289 { 01290 double L = 30.0; // length 01291 double H = 15.0; // height 01292 int ndivx = 30; 01293 int ndivy = 15; 01294 double Lx = L/ndivx; 01295 double Ly = H/ndivy; 01296 01297 // Allocate grid 01298 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed14); 01299 01300 // Allocate material points 01301 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom14, &CbDensity14, &CbAllocMdl14, &CbIniVeloc14, &CbHasTraction14, &CbHasAppDisp14); 01302 01303 // Set callbacks 01304 Prob.pGeom = &CbIsPointInGeom14; 01305 Prob.pB = &CbB14; 01306 Prob.pM = &CbLdM14; 01307 Prob.pVeloc = &CbVel14; 01308 Prob.pStress = &CbStress14; 01309 01310 // Set constants 01311 Prob.t = 0.0; // current time 01312 Prob.tf = 0.7; // final time 01313 Prob.Dt = 0.00001; // time increment 01314 //Prob.tf = 30.0; // final time 01315 //Prob.Dt = 0.001; // time increment 01316 Prob.Dtout = 0.01; // time increment for output 01317 Prob.Outp = 0; // point/particle number to output data 01318 Prob.B = 0.0; // current b-body mass 01319 Prob.M = 0.0; // multiplier for external forces 01320 01321 break; 01322 } 01323 case 15: 01324 { 01325 double L = 10.0; // length 01326 double H = 10.0; // height 01327 int ndivx = 20; 01328 int ndivy = 20; 01329 double Lx = L/ndivx; 01330 double Ly = H/ndivy; 01331 01332 // Allocate grid 01333 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed15); 01334 01335 // Allocate material points 01336 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom15, &CbDensity15, &CbAllocMdl15, &CbIniVeloc15, &CbHasTraction15, &CbHasAppDisp15); 01337 01338 // Set callbacks 01339 Prob.pGeom = &CbIsPointInGeom15; 01340 Prob.pB = &CbB15; 01341 Prob.pM = &CbLdM15; 01342 Prob.pVeloc = &CbVel15; 01343 Prob.pStress = &CbStress15; 01344 01345 // Set constants 01346 Prob.t = 0.0; // current time 01347 Prob.tf = 5.0; // final time 01348 Prob.Dt = 0.0005; // time increment 01349 Prob.Dtout = 0.01; // time increment for output 01350 Prob.Outp = 0; // point/particle number to output data 01351 Prob.B = 0.0; // current b-body mass 01352 Prob.M = 0.0; // multiplier for external forces 01353 01354 break; 01355 } 01356 case 16: 01357 { 01358 double L = 2.0; // length 01359 double H = 10.0; // height 01360 int ndivx = 4; 01361 int ndivy = 20; 01362 double Lx = L/ndivx; 01363 double Ly = H/ndivy; 01364 01365 // Allocate grid 01366 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed16); 01367 01368 // Allocate material points 01369 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom16, &CbDensity16, &CbAllocMdl16, &CbIniVeloc16, &CbHasTraction16, &CbHasAppDisp16); 01370 01371 // Set callbacks 01372 Prob.pGeom = &CbIsPointInGeom16; 01373 Prob.pB = &CbB16; 01374 Prob.pM = &CbLdM16; 01375 Prob.pVeloc = &CbVel16; 01376 Prob.pStress = &CbStress16; 01377 01378 // Set constants 01379 Prob.t = 0.0; // current time 01380 Prob.tf = 0.25; // final time 01381 Prob.Dt = 0.00006; // time increment 01382 Prob.Dtout = 0.001; // time increment for output 01383 Prob.Outp = 0; // point/particle number to output data 01384 Prob.B = 0.0; // current b-body mass 01385 Prob.M = 0.0; // multiplier for external forces 01386 01387 break; 01388 } 01389 case 17: 01390 { 01391 double L = 2.0; // length 01392 double H = 10.0; // height 01393 int ndivx = 4; 01394 int ndivy = 20; 01395 double Lx = L/ndivx; 01396 double Ly = H/ndivy; 01397 01398 // Allocate grid 01399 Prob.G = new Grid2D (/*xMin*/-Lx, /*yMin*/-Ly, /*nRow*/ndivy+3, /*nCol*/ndivx+3, /*Lx*/Lx, /*Ly*/Ly, &CbIsFixed17); 01400 01401 // Allocate material points 01402 Prob.P = new MPoints2D (Prob.NPcell, Prob.G, &CbIsPointInGeom17, &CbDensity17, &CbAllocMdl17, &CbIniVeloc17, &CbHasTraction17, &CbHasAppDisp17); 01403 01404 // Set callbacks 01405 Prob.NPcell = 9; 01406 Prob.pGeom = &CbIsPointInGeom17; 01407 Prob.pB = &CbB17; 01408 Prob.pM = &CbLdM17; 01409 Prob.pVeloc = &CbVel17; 01410 Prob.pStress = &CbStress17; 01411 01412 // Set constants 01413 Prob.t = 0.0; // current time 01414 Prob.tf = 0.25; // final time 01415 Prob.Dt = 0.00006; // time increment 01416 Prob.Dtout = 0.001; // time increment for output 01417 Prob.Outp = 0; // point/particle number to output data 01418 Prob.B = 0.0; // current b-body mass 01419 Prob.M = 0.0; // multiplier for external forces 01420 01421 break; 01422 } 01423 default: { throw new Fatal("Prob # %d is not available.",Prob.Id); } 01424 } 01425 ReSetOutputArrays (Prob); 01426 } 01427 01428 }; // namespace MPM 01429 01430 #endif // MPM_PROBLEMS2D_H
1.7.6.1
