Hi Hugo,
On 15/06/2021 11:28, hugo.levy@gmail.com wrote:
Hi Robert,
I have been using Sfepy for the past two months now, and the more I use it the more I like it! Thank you for developing this great FEM (& more) software.
Thank you, happy to hear that!
I spent some time trying the various functionalities sfepy offers and recently discovered the IGA mode, which I managed to use interactively. After some testing, it turns out isogeometric analysis might be well suited for my applications, where geometry plays a crucial role. So far, my IGA simulations outperform FEM ones, which really sparked my interest for this type of analysis. Having no background in IGA, I read references [1] and [2] to grasp the mathematical concepts at stake, and reference [3] for the Bézier extraction method.
That being said, user's guide mentions that surface terms in IGA are not implemented yet in Sfepy, and I was wondering why?
For instance, I would like to compute integrals over (D-1)-dimensional spaces (where D denotes the problem spatial dimension), but as expected, running the line:
pb.evaluate('ev_surface_integrate.3.Gamma(u)', mode='eval')
results with
ValueError: unknown dof connectivity type! (surface)
This error is raised in sfepy\discrete\iga\fields.py function setup_extra_data.
Of course, I could get around this missing functionality by sampling the scalar field of interest over the desired surface/curve manually (i.e. with field.evaluate_at) and use scipy.integrate. But that would mean giving up on the 'exact' geometry description, which I am not keen on.
Is there a particular reason for not implementing those surface terms in IGA, like a technical difficulty that I am not yet aware of? Otherwise, could you give me insights into where to start if I were to implement my own quadrature rule for computing surface integrals in IGA mode? I haven't figured out in which module the evaluation of Bernstein basis functions at Gauss quadrature points occurs. I know juggling with NURBS, Bézier meshes, etc. doesn't seem particularly easy. Yet I'm willing to give it a try.
There are no unsurmountable technical difficulties - I just quickly implemented the required minimum for our application [4,5], where the domain was cube and no boundary integrals. :) The same reason is behind the single NURBS patch limitation.
If you want to give it a shot, the IGA-related files files are in sfepy/discrete/iga/. For the basis evaluation, there is the initial pure Python implementation in iga.py and more recent Cython implementation in sfepy/discrete/iga/extmods/igac.pyx.
Random notes
In the git version of sfepy, there is a new approximation based on Berstein basis (the same as used in IGA with the Bézier extraction) - it does not allow the smooth boundary as IGA, but the basis is positive, which is important in some applications, and addition of this basis lead to some refactoring, that makes easier the following bullet.
IGA actually evaluates the reference mapping and basis on the boundary when setting the Dirichlet boundary conditions using a l2 projection. The new code is in IGField.get_surface_basis(). I have not thought about it in depth but it should be possible to use something similar for the surface integrals you need.
Supporting several NURBS patches might be technically more involved than adding the surface integrals, because questions regarding mesh format, patch interconnections etc. would have to be sorted out.
There are packages usable from Python that have those capabilities out of the box...
So if you want to give it a try, have a look at the files I mentioned, and do not hesitate to ask more. It would be great to advance the code in this direction.
Cheers, r.
[4] R. Cimrman, M. Novák, R. Kolman, M. Tůma, and J. Vackář, Isogeometric Analysis in Electronic Structure Calculations, Mathematics and Computers in Simulation 145, 125 (2018). [5] R. Cimrman, M. Novák, R. Kolman, M. Tůma, J. Plešek, and J. Vackář, Convergence Study of Isogeometric Analysis Based on Bézier Extraction in Electronic Structure Calculations, Applied Mathematics and Computation 319, 138 (2018).
Thank you very much for your help. Best regards,
Hugo L.
[1] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, Elsevier, 2005, 194 (39-41), pp.4135-4195. 10.1016/j.cma.2004.10.008. hal-01513346
[2] Cimrman, Robert. (2014). Enhancing SfePy with Isogeometric Analysis.
[3] Michael J. Borden, Michael A. Scott, John A. Evans, Thomas J. R. Hughes: Isogeometric finite element data structures based on Bezier extraction of NURBS, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas, March 2010.