3
$\begingroup$

I need to implement a structural analysis code, and I turn to you for advice. My needs are simple:

  1. Library must be integrable in a C++ code
  2. I want to express a weak form, without manual intervention

The second point leaves libraries like deal.II out of my scope. I've tried FEniCS, but it is a nightmare to make it work with Xcode, I am constantly finding discrepancies between console and Xcode, and I need to debug my code!

Any alternatives?

Thanks!

Improve this question
$\endgroup$
4
  • 5
    $\begingroup$ You could ask in the FEniCS support forum for help with Xcode. There is also feelpp.org which offers weak variational forms via TMP. However I have not tried the latter and cannot say whether it is mature enough for everyday use. Another possibility would be to use an arbitrary suitable framework and generate the assemblers using FFC. I did this for two tetrahedral codes and it worked quite well. $\endgroup$
    – Christian Waluga
    Feb 13, 2016 at 17:08
  • $\begingroup$ Thanks @ChristianWaluga, can you tell me more about this last hypothesis with FFC? $\endgroup$
    – senseiwa
    Feb 14, 2016 at 9:48
  • 2
    $\begingroup$ FFC takes a symbolic representation of your weak form and generates all relevant assembly code. Given a suitable finite element code it takes only little effort to plug in the generated local assemblers. This is exactly what FEniCS/Dolfin does. So in case you don't get FEniCS to work within Xcode, you can as well use any other simplex-mesh based framework and still benefit from the elegant FEniCS way of defining your problem. But Wolfgang is right: writing own assemblers requires no magic if you know what you want. It can only be quite error-prone and is not very comfortable for prototyping. $\endgroup$
    – Christian Waluga
    Feb 14, 2016 at 20:17
  • $\begingroup$ Without a more detailed information about the physical nature of your problem, it is difficult to provide a specific advice. If you are looking for an accurate C/C++ code, implementing FEA of shells, the following software (with all the theory described in the applied book) might be useful: members.ozemail.com.au/~comecau/quad_shell.htm $\endgroup$
    – SparseSolverCodes
    Apr 16, 2016 at 11:49

2 Answers 2

Reset to default
5
$\begingroup$

You may not be able to express the weak form in deal.II as a mathematical formula, but you come pretty close. For elasticity, the bilinear form reads $$ a({\mathbf \varphi}_i, {\mathbf \varphi}_j) = \left( \lambda \nabla\cdot {\mathbf \varphi}_i, \nabla\cdot {\mathbf \varphi}_j \right)_\Omega + \left( \mu \nabla\mathbf \varphi_i, \nabla\mathbf \varphi_j \right)_\Omega, + \left( \mu \nabla\mathbf \varphi_i, \nabla\mathbf \varphi_j^T \right)_\Omega, $$ and the corresponding code you would write is as follows:

  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
    for (unsigned int i=0; i<dofs_per_cell; ++i)
      {
        const Tensor<2,dim> phi_i_grad
          = fe_values[displacements].gradient (i,q_point);
        const double phi_i_div
          = fe_values[displacements].divergence (i,q_point);

        for (unsigned int j=0; j<dofs_per_cell; ++j)
          {
            const Tensor<2,dim> phi_j_grad
              = fe_values[displacements].gradient (j,q_point);
            const double phi_j_div
              = fe_values[displacements].divergence (j,q_point);

            cell_matrix(i,j)
              +=  (lambda_values[q_point] *
                   phi_i_div * phi_j_div
                   +
                   mu_values[q_point] *
                   double_contract(phi_i_grad, phi_j_grad)
                   +
                   mu_values[q_point] *
                   double_contract(phi_i_grad, transpose(phi_j_grad))
                  ) *
                  fe_values.JxW(q_point);
          }

The last couple of lines assembling the cell matrix are, in effect, a 1:1 transcription of the formula above.

More details here: https://dealii.org/developer/doxygen/deal.II/group__vector__valued.html

Improve this answer
$\endgroup$
3
$\begingroup$

I myself have not had much experience with computational science. Stil relativity new. However, with the experience that I have had, I will try to answer with the best of my ability.

I would have to agree with the other users here. Try to get FENics support to help troubleshoot your issues with XCODE.

If that does not work, Deal.ii is not that bad. It takes care of much of the hard work of programming everything. With its extensive use of tutorials and documentation, I would not dismiss it quite yet. Sure, it takes a little bit of effort to get up and running, but comparatively speaking, I think that it is one of the best tools there.

As for converting to the weak for and estimating via quadrature, that may take a bit of research.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.