|
00001 // Copyright (C) 2003 Fredrik Bengzon and Johan Jansson. 00002 // Licensed under the GNU GPL Version 2. 00003 00004 #ifndef __ELASTICITY_UPDATED2_H 00005 #define __ELASTICITY_UPDATED2_H 00006 00007 #include <dolfin/PDE.h> 00008 00009 00010 namespace dolfin { 00011 00012 class ElasticityUpdated : public PDE { 00013 public: 00014 00015 ElasticityUpdated(Function::Vector& source, 00016 Function::Vector& wprevious) : PDE(3, 3), f(3), wp(3) 00017 { 00018 add(f, source); 00019 add(wp, wprevious); 00020 } 00021 00022 real lhs(ShapeFunction::Vector& u, ShapeFunction::Vector& v) 00023 { 00024 return 00025 (u(0) * v(0) + u(1) * v(1) + u(2) * v(2)) * dx; 00026 } 00027 00028 real rhs(ShapeFunction::Vector& v) 00029 { 00030 const Point &mp = cell_->midpoint(); 00031 00032 // Debug 00033 00034 /* 00035 // Derivatives 00036 ElementFunction wx0 = wp(0).ddx(); 00037 ElementFunction wy0 = wp(0).ddy(); 00038 ElementFunction wz0 = wp(0).ddz(); 00039 00040 ElementFunction wx1 = wp(1).ddx(); 00041 ElementFunction wy1 = wp(1).ddy(); 00042 ElementFunction wz1 = wp(1).ddz(); 00043 00044 ElementFunction wx2 = wp(2).ddx(); 00045 ElementFunction wy2 = wp(2).ddy(); 00046 ElementFunction wz2 = wp(2).ddz(); 00047 00048 std::cout << "wp(0).ddx(): " << wx0(mp) << std::endl; 00049 std::cout << "wp(0).ddy(): " << wy0(cell_->coord(0)) << std::endl; 00050 std::cout << "wp(0).ddz(): " << wz0(cell_->coord(0)) << std::endl; 00051 00052 std::cout << "wp(1).ddx(): " << wx1(cell_->coord(0)) << std::endl; 00053 std::cout << "wp(1).ddy(): " << wy1(cell_->coord(0)) << std::endl; 00054 std::cout << "wp(1).ddz(): " << wz1(cell_->coord(0)) << std::endl; 00055 00056 std::cout << "wp(2).ddx(): " << wx2(cell_->coord(0)) << std::endl; 00057 std::cout << "wp(2).ddy(): " << wy2(cell_->coord(0)) << std::endl; 00058 std::cout << "wp(2).ddz(): " << wz2(cell_->coord(0)) << std::endl; 00059 */ 00060 00061 // Material parameters 00062 00063 real E = 5.0; 00064 real nu = 0.3; 00065 00066 real lambda = E * nu / ((1 + nu) * (1 - 2 * nu)); 00067 real mu = E / (2 * (1 + nu)); 00068 00069 //real lambda = 20.0; 00070 //real mu = 5.0; 00071 00072 // Compute sigma 00073 00074 static Matrix sigma(3, 3), epsilon(3, 3); 00075 00076 if(!computedsigma[cell_->id()]) 00077 { 00078 00079 static ElementFunction 00080 epsilon00, epsilon01, epsilon02, 00081 epsilon10, epsilon11, epsilon12, 00082 epsilon20, epsilon21, epsilon22; 00083 00084 epsilon(0, 0) = (wp(0).ddx() + wp(0).ddx())(mp); 00085 epsilon(0, 1) = (wp(0).ddy() + wp(1).ddx())(mp); 00086 epsilon(0, 2) = (wp(0).ddz() + wp(2).ddx())(mp); 00087 00088 epsilon(1, 0) = (wp(1).ddx() + wp(0).ddy())(mp); 00089 epsilon(1, 1) = (wp(1).ddy() + wp(1).ddy())(mp); 00090 epsilon(1, 2) = (wp(1).ddz() + wp(2).ddy())(mp); 00091 00092 epsilon(2, 0) = (wp(2).ddx() + wp(0).ddz())(mp); 00093 epsilon(2, 1) = (wp(2).ddy() + wp(1).ddz())(mp); 00094 epsilon(2, 2) = (wp(2).ddz() + wp(2).ddz())(mp); 00095 00096 00097 cout << "epsilon on cell " << cell_->id() << ": " << endl; 00098 epsilon.show(); 00099 00100 sigma(0, 0) = 00101 (lambda + 2 * mu) * epsilon(0, 0) + 00102 lambda * epsilon(1, 1) + lambda * epsilon(2, 2); 00103 sigma(1, 1) = 00104 (lambda + 2 * mu) * epsilon(1, 1) + 00105 lambda * epsilon(0, 0) + lambda * epsilon(2, 2); 00106 sigma(2, 2) = 00107 (lambda + 2 * mu) * epsilon(2, 2) + 00108 lambda * epsilon(0, 0) + lambda * epsilon(1, 1); 00109 sigma(1, 2) = mu * epsilon(1, 2); 00110 sigma(2, 1) = sigma(1, 2); 00111 sigma(2, 0) = mu * epsilon(2, 0); 00112 sigma(0, 2) = sigma(2, 0); 00113 sigma(0, 1) = mu * epsilon(0, 1); 00114 sigma(1, 0) = sigma(0, 1); 00115 00116 //cout << "sigma on cell " << cell_->id() << ": " << endl; 00117 //sigma.show(); 00118 00119 Matrix *sigma0, *sigma1; 00120 00121 sigma0 = sigma0array[cell_->id()]; 00122 sigma1 = sigma1array[cell_->id()]; 00123 00124 sigma *= k; 00125 sigma += *sigma0; 00126 *sigma1 = sigma; 00127 00128 cout << "sigma0 on cell " << cell_->id() << ": " << endl; 00129 sigma0->show(); 00130 cout << "sigma1 on cell " << cell_->id() << ": " << endl; 00131 sigma1->show(); 00132 00133 computedsigma[cell_->id()] = true; 00134 } 00135 else 00136 { 00137 Matrix *sigma1; 00138 sigma1 = sigma1array[cell_->id()]; 00139 00140 sigma = *sigma1; 00141 } 00142 00143 00144 return ((wp(0) * v(0) + wp(1) * v(1) + wp(2) * v(2)) + 00145 k * (f(0) * v(0) + f(1) * v(1) + f(2) * v(2)) - 00146 k * (sigma(0, 0) * v(0).ddx() + 00147 sigma(0, 1) * v(0).ddy() + 00148 sigma(0, 2) * v(0).ddz() + 00149 sigma(1, 0) * v(1).ddx() + 00150 sigma(1, 1) * v(1).ddy() + 00151 sigma(1, 2) * v(1).ddz() + 00152 sigma(2, 0) * v(2).ddx() + 00153 sigma(2, 1) * v(2).ddy() + 00154 sigma(2, 2) * v(2).ddz())) * dx; 00155 00156 //return (wp(0) * v(0) + wp(1) * v(1) + wp(2) * v(2)) * dx; 00157 } 00158 00159 NewArray<Matrix *> sigma0array, sigma1array; 00160 NewArray<bool> computedsigma; 00161 00162 private: 00163 ElementFunction::Vector f; // Source term 00164 ElementFunction::Vector wp; // Velocity value at left end-point 00165 00166 }; 00167 00168 } 00169 00170 #endif 00171 00172 00173 00174
Documentation automatically generated with Doxygen on 10 Sep 2004.