Revision 191cd73559931306f457a6ce9effe70fac906cae authored by Christophe Prud'homme on 19 November 2022, 16:17:50 UTC, committed by Christophe Prud'homme on 19 November 2022, 16:17:50 UTC
remove taskflow trail and prepare for sepcx add support for specx #2009
1 parent 4c07619
qs_elasticity_pure_traction.cpp
/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:set fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4*/
#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feeldiscr/pchv.hpp>
#include <feel/feeldiscr/product.hpp>
#include <feel/feelvf/blockforms.hpp>
#include <feel/feelvf/vf.hpp>
namespace Feel
{
template <typename SpaceType>
NullSpace<double> qsNullSpace( SpaceType const& space, mpl::int_<2> /**/ )
{
auto mode1 = space->element( oneX() );
auto mode2 = space->element( oneY() );
auto mode3 = space->element( vec(Py(),-Px()) );
NullSpace<double> userNullSpace( { mode1,mode2,mode3 } );
return userNullSpace;
}
template <typename SpaceType>
NullSpace<double> qsNullSpace( SpaceType const& space, mpl::int_<3> /**/ )
{
auto mode1 = space->element( oneX() );
auto mode2 = space->element( oneY() );
auto mode3 = space->element( oneZ() );
auto mode4 = space->element( vec(Py(),-Px(),cst(0.)) );
auto mode5 = space->element( vec(-Pz(),cst(0.),Px()) );
auto mode6 = space->element( vec(cst(0.),Pz(),-Py()) );
NullSpace<double> userNullSpace( { mode1,mode2,mode3,mode4,mode5,mode6 } );
return userNullSpace;
}
}
int main(int argc, char**argv )
{
using namespace Feel;
try
{
po::options_description laplacianoptions( "Elasticity options" );
laplacianoptions.add_options()
( "domain", po::value<std::string>()->default_value( "" ), "name of the subdomain where rigid body are applied, if empty then it is the full domain" )
( "E", po::value<double>()->default_value( 1.0e6 ), "Young modulus" )
( "nu", po::value<double>()->default_value( 0.3 ), "Poisson ratio" )
( "no-solve", po::value<bool>()->default_value( false ), "No solve" )
( "weakdir", po::value<bool>()->default_value( false ), "use weak dirichlet" )
( "gamma", po::value<double>()->default_value( 100 ), "penalisation term" )
( "nullspace", po::value<bool>()->default_value( false ), "add null space" )
;
Environment env( _argc=argc, _argv=argv,
_desc=laplacianoptions,
_about=about(_name="qs_elasticity",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
tic();
auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
toc("loadMesh");
tic();
auto Vh = Pchv<1>( mesh );
auto V0hv = Pchv<0>( mesh );
#if FEELPP_DIM == 2
auto V0h = Pch<0>( mesh );
#else
auto V0h = Pchv<0>( mesh );
#endif
auto Xh = product( Vh, V0hv, V0h );
toc("Vh");
auto u = Vh->element("u");
auto v = Vh->element("v");
auto mv = V0hv->element("mv");
auto nv = V0hv->element("nv");
auto m = V0h->element("m");
auto n = V0h->element("n");
auto nu = doption(_name="nu");
auto E = doption(_name="E");
auto lambda = E*nu/( (1+nu)*(1-2*nu) );
auto mu = E/(2*(1+nu));
auto deft = sym(gradt(u));
auto def = sym(grad(u));
auto Id = eye<FEELPP_DIM,FEELPP_DIM>();
auto sigmat = lambda*trace(deft)*Id + 2*mu*deft;
auto sigma = lambda*trace(def)*Id + 2*mu*def;
auto f = expr<FEELPP_DIM,1>( soption(_name="functions.f"), "f" );
auto g = expr<FEELPP_DIM,1>( soption(_name="functions.g"), "g" );
auto tau = expr<FEELPP_DIM,1>( soption(_name="functions.tau"), "tau" );
#if FEELPP_DIM == 2
auto omega = expr( soption(_name="functions.omega"), "omega" );
#else
auto omega = expr<3,1>( soption(_name="functions.omega"), "omega" );
#endif
auto rigid_body_domain_name = soption(_name="domain");
tic();
auto l = blockform1( Xh, solve::strategy::monolithic, backend() );
l(0_c) = integrate(_range=elements(mesh),
_expr=inner(f,id(v)));
l(0_c) += integrate(_range=markedfaces(mesh,"Neumann"),
_expr=inner(g,id(v)));
auto rigid_body_domain = (rigid_body_domain_name.empty())?elements(mesh):markedelements(mesh, rigid_body_domain_name );
std::cout << "surface rigid domain=" << integrate( _range=rigid_body_domain, _expr=cst(1.0) ).evaluate() << std::endl;
// translation
l(1_c) += integrate(_range=elements(mesh),
_expr=inner(tau,id(mv)));
// rotation
l(2_c) += integrate(_range=rigid_body_domain,
_expr=inner(omega,id(m)));
toc("l");
tic();
auto a = blockform2( Xh, solve::strategy::monolithic, backend() );
a(0_c, 0_c) = integrate(_range=elements(mesh),
_expr=inner( sigmat, grad(v) ) );
a(0_c,1_c) += integrate(_range=elements(mesh),
_expr=trans(idt(mv))*id(v));
a(1_c,0_c) += integrate(_range=elements(mesh),
_expr=trans(id(mv))*idt(v));
#if FEELPP_DIM == 2
a(0_c,2_c) += integrate(_range=rigid_body_domain,
_expr=trans(idt(m))*curlx(v));
a(2_c,0_c) += integrate(_range=rigid_body_domain,
_expr=trans(id(m))*curlxt(v));
#else
a(0_c,2_c) += integrate(_range=rigid_body_domain,
_expr=trans(idt(m))*curl(v));
a(2_c,0_c) += integrate(_range=rigid_body_domain,
_expr=trans(id(m))*curlt(v));
#endif
toc("a");
//! solve the linear system, find u s.t. a(u,v)=l(v) for all v
if ( !boption( "no-solve" ) )
{
tic();
std::shared_ptr<NullSpace<double> > myNullSpace( new NullSpace<double>(backend(),qsNullSpace(Vh,mpl::int_<FEELPP_DIM>())) );
backend()->attachNearNullSpace( myNullSpace );
if ( boption(_name="nullspace") )
backend()->attachNearNullSpace( myNullSpace );
auto U=Xh.element();
a.solve(_rhs=l,_solution=U);
u=U(0_c);
//m=U(1_c);
toc("a.solve");
}
tic();
auto e = exporter( _mesh=mesh );
e->addRegions();
e->add( "u", u );
e->save();
toc("Exporter");
return 0;
}
catch( ... )
{
handleExceptions();
}
return 1;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...