On 2 Dec 2016, at 15:01, Bruno Turcksin <[email protected]> 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.
