Commit 007ab19c authored by Alexandre Halbach's avatar Alexandre Halbach

Add fully anisotropic elasticity

parent 83f7b412
......@@ -32,7 +32,7 @@ void sparselizard(void)
formulation elasticity;
// The linear elasticity formulation is classical and thus predefined:
elasticity += integral(vol, predefinedelasticity(u, E, nu));
elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), E, nu));
// Add the inertia terms:
elasticity += integral(vol, -rho*dtdt(dof(u))*tf(u));
......
......@@ -32,7 +32,7 @@ void sparselizard(void)
formulation elasticity;
// The linear elasticity formulation is classical and thus predefined:
elasticity += integral(vol, predefinedelasticity(u, E, nu));
elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), E, nu));
// Add a volumic force in the -z direction:
elasticity += integral(vol, array1x3(0,0,-10)*tf(u));
......
......@@ -126,7 +126,7 @@ void sparselizard(void)
electrostatics += integral(electricdomain, 20, epsilon*transpose(invjac*grad(dof(v)))*invjac*grad(tf(v)) * detjac );
// The linear elasticity formulation is classical and thus predefined:
elastoacoustic += integral(solid, predefinedelasticity(u, E, nu, "planestrain"));
elastoacoustic += integral(solid, predefinedelasticity(dof(u), tf(u), E, nu, "planestrain"));
// Add the inertia terms:
elastoacoustic += integral(solid, -rho*dtdt(dof(u))*tf(u));
......
......@@ -79,7 +79,7 @@ void sparselizard(void)
electrostatics += integral(electricdomain, umesh, epsilon*grad(dof(v))*grad(tf(v)));
// The linear elasticity formulation is classical and thus predefined:
elasticity += integral(solid, predefinedelasticity(u, E, nu, "planestrain"));
elasticity += integral(solid, predefinedelasticity(dof(u), tf(u), E, nu, "planestrain"));
// Electrostatic forces, computed on the elements of the whole electric domain
// but with mechanical deflection test functions tf(u) only on solid.
......
......@@ -341,73 +341,57 @@ expression mathop::m2d(expression input)
return expression(3,1,{dx(compx(input)), dy(compy(input)), dy(compx(input)) + dx(compy(input))});
}
expression mathop::m3dn(expression input)
expression mathop::m3d(expression input)
{
if (input.countrows() != 3 || input.countcolumns() != 1)
{
std::cout << "Error in 'mathop' namespace: can only compute m3dn of a 3x1 column vector" << std::endl;
std::cout << "Error in 'mathop' namespace: can only compute m3d of a 3x1 column vector" << std::endl;
abort();
}
return expression(3,1,{dx(compx(input)), dy(compy(input)), dz(compz(input))});
return expression(6,1,{dx(compx(input)), dy(compy(input)), dz(compz(input)), dz(compy(input)) + dy(compz(input)), dx(compz(input)) + dz(compx(input)), dx(compy(input)) + dy(compx(input))});
}
expression mathop::m3ds(expression input)
////////// PREDEFINED FORMULATIONS
expression mathop::predefinedelasticity(expression dofu, expression tfu, expression E, expression nu, std::string myoption)
{
if (input.countrows() != 3 || input.countcolumns() != 1)
{
std::cout << "Error in 'mathop' namespace: can only compute m3ds of a 3x1 column vector" << std::endl;
abort();
}
return expression(3,1,{dx(compy(input)) + dy(compx(input)), dz(compy(input)) + dy(compz(input)), dx(compz(input)) + dz(compx(input))});
expression hookesmatrix(6,6, {1-nu,nu,nu,0,0,0, nu,1-nu,nu,0,0,0, nu,nu,1-nu,0,0,0, 0,0,0,0.5*(1-2*nu),0,0, 0,0,0,0,0.5*(1-2*nu),0, 0,0,0,0,0,0.5*(1-2*nu)});
hookesmatrix = E/(1+nu)/(1-2*nu) * hookesmatrix;
return predefinedelasticity(dofu, tfu, hookesmatrix, myoption);
}
////////// PREDEFINED FORMULATIONS
expression mathop::predefinedelasticity(expression u, expression E, expression nu, std::string myoption)
expression mathop::predefinedelasticity(expression dofu, expression tfu, expression hookesmatrix, std::string myoption)
{
if (u.countrows() == 1)
if (dofu.countrows() != tfu.countrows() || dofu.countcolumns() != 1 || tfu.countcolumns() != 1 || dofu.countrows() == 1)
{
std::cout << "Error in 'mathop' namespace: 'predefinedelasticity' is undefined in 1D" << std::endl;
std::cout << "Error in 'mathop' namespace: first two arguments in 'predefinedelasticity' must be both either 2x1 or 3x1 vectors" << std::endl;
abort();
}
if (u.countrows() == 2)
if (dofu.countrows() == 2)
{
// There must be an option in 2D:
if (myoption.length() == 0)
{
std::cout << "Error in 'mathop' namespace: for a 2D problem a last string argument is needed in 'predefinedelasticity'" << std::endl;
std::cout << "Available choices are: 'planestrain', 'planestress', 'axisymmetry'" << std::endl;
abort();
}
if (myoption == "planestrain")
return - transpose(E/(1+nu)/(1-2*nu)*array3x3(1-nu,nu,0,nu,1-nu,0,0,0,0.5*(1-2*nu))*m2d(dof(u)))*m2d(tf(u));
if (myoption == "planestress")
return - transpose(E/(1-pow(nu,2))*array3x3(1,nu,0,nu,1,0,0,0,0.5*(1-nu))*m2d(dof(u)))*m2d(tf(u));
if (myoption == "axisymmetry")
{
// The radius in cylindrical coordinates:
field x("x");
expression H = E/(1+nu)/(1-2*nu) * expression(4,4,{1-nu,nu,nu,0, nu,1-nu,nu,0, nu,nu,1-nu,0, 0,0,0,0.5-nu});
expression H = expression(3,3, {
hookesmatrix.getarrayentry(0,0),hookesmatrix.getarrayentry(0,1),hookesmatrix.getarrayentry(0,5),
hookesmatrix.getarrayentry(1,0),hookesmatrix.getarrayentry(1,1),hookesmatrix.getarrayentry(1,5),
hookesmatrix.getarrayentry(5,0),hookesmatrix.getarrayentry(5,1),hookesmatrix.getarrayentry(5,5) });
expression dofexpr = expression(4,1,{dx(compx(dof(u))),dy(compy(dof(u))),1/x*compx(dof(u)),dx(compy(dof(u)))+dy(compx(dof(u)))});
expression tfexpr = expression(4,1, {dx(compx(tf(u))), dy(compy(tf(u))), 1/x*compx(tf(u)), dx(compy(tf(u))) +dy(compx(tf(u)))});
// Also multiply by the coordinate change Jacobian determinant, i.e. by the radius (x here):
return - transpose(H*dofexpr)*tfexpr * x;
return -transpose( H *m2d(dofu) )*m2d(tfu);
}
// If the option is not valid:
std::cout << "Error in 'mathop' namespace: invalid option or no option provided for the 2D problem in 'predefinedelasticity'" << std::endl;
std::cout << "Available choices are: 'planestrain'" << std::endl;
abort();
}
if (u.countrows() == 3)
if (dofu.countrows() == 3)
{
if (myoption.length() > 0)
{
std::cout << "Error in 'mathop' namespace: for a 3D problem the last string argument must be empty in 'predefinedelasticity'" << std::endl;
abort();
}
return - transpose(E/(1+nu)/(1-2*nu)*array3x3(1-nu,nu,nu,nu,1-nu,nu,nu,nu,1-nu)*m3dn(dof(u)))*m3dn(tf(u)) - transpose(E/(1+nu)/(1-2*nu)*array3x3(0.5*(1-2*nu),0,0,0,0.5*(1-2*nu),0,0,0,0.5*(1-2*nu))*m3ds(dof(u)))*m3ds(tf(u));
return -transpose( hookesmatrix*m3d(dofu) ) * m3d(tfu);
}
}
......
......@@ -111,12 +111,15 @@ namespace mathop
////////// PREDEFINED OPERATORS
expression m2d(expression input);
expression m3dn(expression input);
expression m3ds(expression input);
expression m3d(expression input);
////////// PREDEFINED FORMULATIONS
expression predefinedelasticity(expression u, expression Eyoung, expression nupoisson, std::string myoption = "");
// Isotropic linear elasticity:
expression predefinedelasticity(expression dofu, expression tfu, expression Eyoung, expression nupoisson, std::string myoption = "");
// General anisotropic linear elasticity:
expression predefinedelasticity(expression dofu, expression tfu, expression hookesmatrix, std::string myoption = "");
expression predefinedelectrostaticforce(expression gradtfu, expression gradv, expression epsilon);
};
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment