Thanks, this is really interesting!
Franco
On Monday, December 5, 2016 at 11:35:35 AM UTC+1, Luca Heltai wrote:
>
> On 2 Dec 2016, at 15:01, Bruno Turcksin <[email protected]
> <javascript:>> wrote:
> >
> >
> > 2 - Are there facilities of some kind that can help in generating the
> assembly code? In Fenics I just specified the weak formulations.
> > You need to write the code in C++ and loop over cells, quadrature
> points, and basis functions yourself. So this will be more verbose than
> what you do in Fenics.
> >
>
>
> However, there are projects that will help you prototype your application
> (in threadsafe/parallel mode), such as pi-domus
>
> https://github.com/mathLab/pi-DoMUS
>
> Here you can solve a non-linear, time dependent problem, by specifying
> just the energy of your system (similarly to what is done in the fenics).
> All system matrices are computed in parallel, using both MPI and TBB, by
> automatic differentiation.
>
> See, for example, the interface
>
>
> https://github.com/mathLab/pi-DoMUS/blob/master/include/interfaces/neo_hookean_two_fields.h
>
>
> Interface files are the only thing needed to have the problem solve
> non-linear time dependent pdes in pi-domus, e.g.:
>
>
>
> template <int dim, int spacedim,typename LAC>
> template<typename EnergyType,typename ResidualType>
> void
> NeoHookeanTwoFieldsInterface<dim,spacedim,LAC>::energies_and_residuals(const
> typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
> FEValuesCache<dim,spacedim> &fe_cache,
> std::vector<EnergyType> &energies,
> std::vector<std::vector<ResidualType> > &,
> bool compute_only_system_terms) const
> {
> EnergyType alpha = 0;
> this->reinit(alpha, cell, fe_cache);
>
> const FEValuesExtractors::Vector displacement(0);
>
> auto &grad_us = fe_cache.get_gradients("solution", "u", displacement,
> alpha);
> auto &Fs = fe_cache.get_deformation_gradients("solution", "Fu",
> displacement, alpha);
>
> const FEValuesExtractors::Scalar pressure(dim);
> auto &ps = fe_cache.get_values("solution","p", pressure, alpha);
>
> auto &JxW = fe_cache.get_JxW_values();
>
> const unsigned int n_q_points = ps.size();
>
> energies[0] = 0;
> for (unsigned int q=0; q<n_q_points; ++q)
> {
> Tensor <1, dim, double> B;
>
> const EnergyType &p = ps[q];
> const Tensor <2, dim, EnergyType> &F = Fs[q];
> const Tensor<2, dim, EnergyType> C = transpose(F)*F;
> auto &grad_u = grad_us[q];
>
> EnergyType Ic = trace(C);
> EnergyType J = determinant(F);
>
> EnergyType psi = (mu/2.)*(Ic-dim)-p*(J-1.);
> energies[0] += (psi)*JxW[q];
> if (!compute_only_system_terms)
> energies[1] += (scalar_product(grad_u,grad_u) + 0.5*p*p)*JxW[q];
> }
> }
>
>
> with the above 60~ lines of codes you can solve a non-linear
> incompressible elasticity problem using two fields.
>
> L.
>
>
>
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
