Implementation of IGA surface terms
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.
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.
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.
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.
Thank you very much for your swift reply! I will take a closer look at the code, now that I know where the relevant part is located. I will post updates on that discussion thread.
Cheers,
Hugo L.
Hi Robert,
So, I took an in-depth look at the various routines defined in sfepy\discrete\iga\iga.py sfepy\discrete\iga\field.py as well as the Cython implementation in extmods. As you mentioned, I found myself particularly interested in IGField.get_surface_basis() which seems to provide all the input required to compute the surface/curve integral.
Down below is how I slightly adapted the existing code to my needs (comments included):
def get_surface_integral(self, region, U): """ U : 1D numpy array scalar field iga-dofs defined over the whole domain (omega) """ from sfepy.discrete.integrals import Integral import sfepy.discrete.iga as iga # from sfepy.discrete.iga.extmods.igac import eval_mapping_data_in_qp from sfepy.discrete.iga.iga import eval_mapping_data_in_qp # use the python implementation to understand what's going on under the hoods
nurbs = self.nurbs
facets = self._get_facets(region.kind_tdim)
# assert scalar field U has the same length as control points
assert U.shape[0] == nurbs.cps.shape[0], "scalar field and control points do not have the same length"
# Cell and face(cell) ids for each facet.
fis = region.get_facet_indices()
# Integral given by max. NURBS surface degree.
fdegrees = iga.get_surface_degrees(nurbs.degrees)
order = fdegrees.max()
integral = Integral('i', order=2*order)
vals, weights = integral.get_qp(self.domain.gel.surface_facet_name)
# Boundary QP - use tensor product structure.
bvals = iga.create_boundary_qp(vals, region.tdim)
# Compute facet basis, jacobians and physical BQP.
intU = 0.0 # integral value initialization
for ii, (ie, ifa) in enumerate(fis):
qp_coors = bvals[ifa]
bfs, _, dets = eval_mapping_data_in_qp(qp_coors, nurbs.cps,
nurbs.weights,
nurbs.degrees,
nurbs.cs,
nurbs.conn,
nm.array([ie]))
# Facet basis.
fbfs = bfs[..., facets[ifa]][0, :, 0, :]
# Weight Jacobians by quadrature point weights.
dets = nm.abs(dets) * weights[None, :, None, None]
dets = dets[0, :, 0, :]
# Physical BQP.
#---- code modification ----#
U_fcps = U[nurbs.conn[ie, facets[ifa]]] # control points replaced by scalar field
U_qp = nm.dot(fbfs, U_fcps) # compute the value of the scalar field at quadrature points
seg = nm.dot(U_qp, dets) # quadrature rule
intU += seg
return intUand I am calling this IGField method like this (the domain is defined by two concentric circles of radii 0.5 and 1.0):
from sfepy.discrete.iga.domain import IGDomain from sfepy.discrete.common.fields import Field import numpy as np
# Function definitions def get_inner(coors, domain=None): x, y = coors[:, 0], coors[:, 1]
r = np.sqrt(x**2 + y**2)
flag = np.where(( abs(r-0.5) < 1e-10 ))[0]
return flag# Setting a Domain instance domain = IGDomain.from_file('circles.iga')
# Sub-domains omega = domain.create_region('Omega', 'all') Gamma_in = domain.create_region('Gamma_in', 'vertices by get_inner', kind='facet', functions={'get_inner' : get_inner})
# field
field = Field.from_args('fu', np.float64, 'scalar', omega,
approx_order=None, space='H1', poly_space_base='iga')
nurbs = domain.nurbs n_dofs = len(nurbs.cps) U = np.ones(n_dofs) intU = field.get_surface_integral(Gamma_in, U)
Unfortunately, this does not work. Here is what I understand:
- the scalar field evaluation at quadrature points U_qp seems to return the correct result (from my tests).
- for a scalar field constant over Gamma_in (say U is 1.0 on Gamma_in), get_surface_integral() should return the inner circle perimeter, that is pi. Given the previous bullet, I infer that the weights used to define the quadrature rule on a single '1_2' element is not correct.
Variable 'dets' takes into account both the jacobians of the mapping from reference bezier element (i.e. abs( det( d_x / d_xi ) ) if I am not mistaken) AND the weights associated with the quadrature points for '1_2' line elements.
I might be missing something, but haven't figured out what yet... I would be grateful if you could help me with this. I can share the 'circles.iga' file in case (which was generated with the create_from_igakit and write_iga_data routines).
Thanks in advance for your time. Cheers,
Hugo L.
Hi Hugo,
On 17/06/2021 11:37, hugo.levy@gmail.com wrote:
Hi Robert,
So, I took an in-depth look at the various routines defined in sfepy\discrete\iga\iga.py sfepy\discrete\iga\field.py as well as the Cython implementation in extmods. As you mentioned, I found myself particularly interested in IGField.get_surface_basis() which seems to provide all the input required to compute the surface/curve integral.
Ultimately, IGField.setup_extra_data() will have to be properly implemented, and surface IGMapping created in IGField.create_mapping(). But your test is the first step in the right direction.
<snip>
Unfortunately, this does not work. Here is what I understand:
- the scalar field evaluation at quadrature points U_qp seems to return the correct result (from my tests).
- for a scalar field constant over Gamma_in (say U is 1.0 on Gamma_in), get_surface_integral() should return the inner circle perimeter, that is pi. Given the previous bullet, I infer that the weights used to define the quadrature rule on a single '1_2' element is not correct.
Variable 'dets' takes into account both the jacobians of the mapping from reference bezier element (i.e. abs( det( d_x / d_xi ) ) if I am not mistaken) AND the weights associated with the quadrature points for '1_2' line elements.
I might be missing something, but haven't figured out what yet... I would be grateful if you could help me with this. I can share the 'circles.iga' file in case (which was generated with the create_from_igakit and write_iga_data routines).
Yes, please, send me the file so that I can run the code (off-list if needed).
Cheers, r.
Hi Robert,
I have finally found where the error lies in the approach I proposed few days ago:
bfs, _, dets = eval_mapping_data_in_qp(qp_coors, nurbs.cps, nurbs.weights, nurbs.degrees, nurbs.cs, nurbs.conn, nm.array([ie]))
This part is wrong. More specifically, the jacobian determinants are not the correct ones for performing a surface/line integral. Let me try to explain.
** Volume Integration ** When performing the integration over the whole domain Omega, one needs to compute the Jacobian associated with the mapping (in 2D) (\xi, \eta) --> (x, y). The Jacobian matrix is therefore a 2*2 matrix J = [ [dx/dxi dx/deta] , [dy/dxi dy/deta] ]
** Surface/line Inegration ** Here, the mapping is actually different and is of the form sigma --> (\xi(sigma) , \eta(sigma)) --> (x, y) ; where (x, y) belongs to Gamma. In my case, and because Gamma is actually a part of the boundary, its representation in the parameter space corresponds to the left border (\xi=0 , eta). It follows that the Jacobian of the transformation is just J = [dx/deta , dy/deta] (which is simply the second column of the full matrix written above).
With the correct jacobian, I am able to retrieve the correct perimeter for a circle of radius 0.5, that is pi.
Now I thinking on how to generalize a bit more this approach (i.e. not having to identify the boundary in the parameter space 'by hand'). As always, I am open to any suggestion or remark.
Thanks for your help!
Hugo L.
Hi Hugo,
On 21/06/2021 14:18, hugo.levy@gmail.com wrote:
Hi Robert,
I have finally found where the error lies in the approach I proposed few days ago:
bfs, _, dets = eval_mapping_data_in_qp(qp_coors, nurbs.cps, nurbs.weights, nurbs.degrees, nurbs.cs, nurbs.conn, nm.array([ie]))
This part is wrong. More specifically, the jacobian determinants are not the correct ones for performing a surface/line integral. Let me try to explain.
** Volume Integration ** When performing the integration over the whole domain Omega, one needs to compute the Jacobian associated with the mapping (in 2D) (\xi, \eta) --> (x, y). The Jacobian matrix is therefore a 2*2 matrix J = [ [dx/dxi dx/deta] , [dy/dxi dy/deta] ]
** Surface/line Inegration ** Here, the mapping is actually different and is of the form sigma --> (\xi(sigma) , \eta(sigma)) --> (x, y) ; where (x, y) belongs to Gamma. In my case, and because Gamma is actually a part of the boundary, its representation in the parameter space corresponds to the left border (\xi=0 , eta). It follows that the Jacobian of the transformation is just J = [dx/deta , dy/deta] (which is simply the second column of the full matrix written above).
You are right, the differential line element is d = [dx/dsigma dy/dsigma], and the line element size for integration is |d| (|| = l2 norm).
So it seems that also IGField.get_surface_basis() has the area determinant wrong, as it uses the one from the volume mapping, the basis functions seem OK. -> You found a bug that should manifest when setting non-constant Dirichlet boundary conditions on a domain that is not just a rectangle/box aligned with the coordinate axes. (*)
Clearly, the IGA code is a bit rusty, and would benefit from a face lift :)
With the correct jacobian, I am able to retrieve the correct perimeter for a circle of radius 0.5, that is pi.
With the wrong implementation, you get pi/2 instead of pi, right? (That is what I get with the code you sent me.)
Now I thinking on how to generalize a bit more this approach (i.e. not having to identify the boundary in the parameter space 'by hand'). As always, I am open to any suggestion or remark.
Random thoughts:
To make in properly, one would need to implement an analogous function to _s_describe() in sfepy/discrete/common/extmods/refmaps.c into igac.pyx, that computes also the surface normal. The surface integral is always(?) over a boundary side of the NURBS patch, that is always a rectangle in the parametric coordinates. This should allow to select correct volume Jacobian matrix components for the surface element mapping "area element magnitude", based on facet indices. See also get_patch_box_regions(), get_bezier_element_entities() and IGField.get_econn().
One could construct explicitely the 2D in 3D or 1D in 2D NURBS patch using the tensor product structure, but then one must be careful to have compatible orientations and numbering with the corresponding facet of the underlying 3D or 2D patch, respectively. So IMO taking the subset of volume basis that corresponds to a face is an easier way. In FEM, it is what we do now (in distant past the surface basis functions were constructed independently, but that lead to subtle ordering/numbering bugs...)
Currently there is the single NURBS patch limitation in the code, so from IGA point of view, only a single NUBRS element is present - what we see is the auxiliary mesh of Bezier elements for making the Bezier extraction operator (as in the references you mentioned earlier). That is why get_patch_box_regions() in iga.py works. But that probably is not important for the surface integrals code, besides the fact that it would not be possible to integrate over inner domain surfaces (e.g. material boundaries) without the multi-patch support.
I do not have currently time to dive into it too deeply, but I will try to respond in a timely manner, should you decide to go on!
Cheers, r.
(*) It would be nice to verify that this bug is indeed there.
Hi Robert,
You are right, the differential line element is d = [dx/dsigma dy/dsigma], and the line element size for integration is |d| (|| = l2 norm).
Yes, that's clearer put that way. The volume-mapping-jacobian-matrix (let's call it J) can be re-used for computing [dx/dsigma dy/dsigma] = J * [dxi/dsigma deta/dsigma]
The surface integral is always(?) over a boundary side of the NURBS patch, that is always a rectangle in the parametric coordinates. This should allow to select correct volume Jacobian matrix components for the surface element mapping "area element magnitude", based on facet indices.
I believe you are right. In my case for instance, I see that variable ifa=3 which means the line-integration is carried out over the left-side of the reference square element. Thus [dxi/dsigma deta/dsigma] = [dxi/deta deta/deta] = [0 1]. If I want to integrate a scalar field over a random line defined on Omega (i.e. not necessarly a part of the boundary), this property is lost. Fortunately, this is a case of limited interest for me, which means I will never have to actually compute [dxi/dsigma deta/dsigma], but rather know which side of the square is at stake. As you mention, it would be easier with the multi-patch support.
Some thoughts / questions:
For a surface integral in 3D, one needs to introduce a second parameter lambda in addition to sigma. The differential surface element is then s = [dx/dsigma dy/dsigma] cross_product [dx/dlambda dy/dlambda]. The reference element is the unit cube. However, I did not find the numbering of the faces used in Sfepy...
I have not thought of the case where the integral is carried out over a facet / edge belonging to two cells. But again, I believe that case never occurs for a boundary integral.
With the wrong implementation, you get pi/2 instead of pi, right? (That is what I get with the code you sent me.)
With the 'circles_coarse.iga' file that I sent you, yes, I get roughly pi/2. This is actually an odd coincidence as I got completely different results when refining the iga domain with more knots.
One could construct explicitely the 2D in 3D or 1D in 2D NURBS patch using the tensor product structure, but then one must be careful to have compatible orientations and numbering with the corresponding facet of the underlying 3D or 2D patch, respectively.
This is the very first idea I had (kind of), then realized I did not know how to extract the correct dofs. However, I gained a lot of insight into how iga in sfepy works and might give it another try now. To test this concept, I will generate two separate .iga files and see if that works. Thanks for the idea.
You found a bug that should manifest when setting non-constant Dirichlet boundary conditions on a domain that is not just a rectangle/box aligned with the coordinate axes. (*) (*) It would be nice to verify that this bug is indeed there.
I could test that. However I do not see when IGField.get_surface_basis() is called? Not in pb.set_bcs() apparently.
Many thanks for your help. At least I can now do the computation that I originally asked help for, which is great.
Cheers,
Hugo L.
Hi Hugo,
On 22/06/2021 11:00, hugo.levy@gmail.com wrote:
Hi Robert,
You are right, the differential line element is d = [dx/dsigma dy/dsigma], and the line element size for integration is |d| (|| = l2 norm).
Yes, that's clearer put that way. The volume-mapping-jacobian-matrix (let's call it J) can be re-used for computing [dx/dsigma dy/dsigma] = J * [dxi/dsigma deta/dsigma]
The surface integral is always(?) over a boundary side of the NURBS patch, that is always a rectangle in the parametric coordinates. This should allow to select correct volume Jacobian matrix components for the surface element mapping "area element magnitude", based on facet indices.
I believe you are right. In my case for instance, I see that variable ifa=3 which means the line-integration is carried out over the left-side of the reference square element. Thus [dxi/dsigma deta/dsigma] = [dxi/deta deta/deta] = [0 1]. If I want to integrate a scalar field over a random line defined on Omega (i.e. not necessarly a part of the boundary), this property is lost. Fortunately, this is a case of limited interest for me, which means I will never have to actually compute [dxi/dsigma deta/dsigma], but rather know which side of the square is at stake. As you mention, it would be easier with the multi-patch support.
For a random line integral, you could sample values along the curve using FieldVariable.evaluate_at(). (Analogously in 3D.)
Some thoughts / questions:
For a surface integral in 3D, one needs to introduce a second parameter lambda in addition to sigma. The differential surface element is then s = [dx/dsigma dy/dsigma] cross_product [dx/dlambda dy/dlambda]. The reference element is the unit cube. However, I did not find the numbering of the faces used in Sfepy...
See sfepy/discrete/fem/geometry_element.py, the geometry_data dict. This information is accessible through IDdomain.gel. The ifa above relates to the row index into the faces array.
I have not thought of the case where the integral is carried out over a facet / edge belonging to two cells. But again, I believe that case never occurs for a boundary integral.
Usually it is not desirable to integrate over the common (internal) facet twice, that is why SfePy allows to define a parent region of a region, and then it is ensured that the cell/facet combination with normal going outside(IIRC) the parent region is used [1].
[1] https://sfepy.org/doc-devel/examples/linear_elasticity-linear_elastic_tracti...
With the wrong implementation, you get pi/2 instead of pi, right? (That is what I get with the code you sent me.)
With the 'circles_coarse.iga' file that I sent you, yes, I get roughly pi/2. This is actually an odd coincidence as I got completely different results when refining the iga domain with more knots.
Indeed, this is quite misleading :)
One could construct explicitely the 2D in 3D or 1D in 2D NURBS patch using the tensor product structure, but then one must be careful to have compatible orientations and numbering with the corresponding facet of the underlying 3D or 2D patch, respectively.
This is the very first idea I had (kind of), then realized I did not know how to extract the correct dofs. However, I gained a lot of insight into how iga in sfepy works and might give it another try now. To test this concept, I will generate two separate .iga files and see if that works. Thanks for the idea.
You found a bug that should manifest when setting non-constant Dirichlet boundary conditions on a domain that is not just a rectangle/box aligned with the coordinate axes. (*) (*) It would be nice to verify that this bug is indeed there.
I could test that. However I do not see when IGField.get_surface_basis() is called? Not in pb.set_bcs() apparently.
pb.set_bcs() only deals with the BCs definitions. The actual calculation of BCs' values occurs in Equations.time_update() -> Variables.equation_mapping() -> EquationMap.map_equations() (sfepy/discrete/common/dof_info.py), where Field.set_dofs() (sfepy/discrete/common/fields.py) is called on the user-provided data. There: from sfepy.discrete.projections import project_to_facets, and finally project_to_facets() calls field.get_surface_basis(region). It's a bit convoluted :).
BTW. the test could be just a problem description file and a non-box-like domain (e.g. your circle), where you set ebcs by a function, and in a postprocess() function call variable.evaluate_at() on the boundary to verify the values provided.
Many thanks for your help. At least I can now do the computation that I originally asked help for, which is great.
Hth!
r.
Hi Robert,
See sfepy/discrete/fem/geometry_element.py, the geometry_data dict. This information is accessible through IDdomain.gel. The ifa above relates to the row index into the faces array.
Thanks! Indeed I was looking for info in the discrete/iga directory exclusively. Now the surface integral in 3D also works flawlessly .
pb.set_bcs() only deals with the BCs definitions. The actual calculation of BCs' values occurs in Equations.time_update() -> Variables.equation_mapping() -> EquationMap.map_equations() (sfepy/discrete/common/dof_info.py), where Field.set_dofs() (sfepy/discrete/common/fields.py) is called on the user-provided data. There: from sfepy.discrete.projections import project_to_facets, and finally project_to_facets() calls field.get_surface_basis(region). It's a bit convoluted :).
BTW. the test could be just a problem description file and a non-box-like domain (e.g. your circle), where you set ebcs by a function, and in a postprocess() function call variable.evaluate_at() on the boundary to verify the values provided.
Okay, that's quite intricate as you say. I did the test you suggested by setting the dirichlet condition on the inner circle to the x-coordinate. I then solved a laplace pde and called field.evaluate_at() at various locations on the inner circle. No bugs observed here, the values I got were the expected ones. Again, I can send you the problem description file I used in case you want to verify by yourself.
Anyway, thank you for your help and patience. I hope this discussion thread will be helpful to others eventually. I will definitely continue to work with sfepy for my IGA-related applications!
Cheers,
Hugo
Hi Hugo,
On 23/06/2021 16:27, hugo.levy@gmail.com wrote:
Hi Robert,
See sfepy/discrete/fem/geometry_element.py, the geometry_data dict. This information is accessible through IDdomain.gel. The ifa above relates to the row index into the faces array.
Thanks! Indeed I was looking for info in the discrete/iga directory exclusively. Now the surface integral in 3D also works flawlessly .
pb.set_bcs() only deals with the BCs definitions. The actual calculation of BCs' values occurs in Equations.time_update() -> Variables.equation_mapping() -> EquationMap.map_equations() (sfepy/discrete/common/dof_info.py), where Field.set_dofs() (sfepy/discrete/common/fields.py) is called on the user-provided data. There: from sfepy.discrete.projections import project_to_facets, and finally project_to_facets() calls field.get_surface_basis(region). It's a bit convoluted :).
BTW. the test could be just a problem description file and a non-box-like domain (e.g. your circle), where you set ebcs by a function, and in a postprocess() function call variable.evaluate_at() on the boundary to verify the values provided.
Okay, that's quite intricate as you say. I did the test you suggested by setting the dirichlet condition on the inner circle to the x-coordinate. I then solved a laplace pde and called field.evaluate_at() at various locations on the inner circle. No bugs observed here, the values I got were the expected ones. Again, I can send you the problem description file I used in case you want to verify by yourself.
Yes, that would be great. It would be helpful to have for the eventual proper implementation of surface integrals.
Anyway, thank you for your help and patience. I hope this discussion thread will be helpful to others eventually. I will definitely continue to work with sfepy for my IGA-related applications!
You are welcome! If you have a nice IGA-related example, that could replace or complement the current ones, you may consider making a pull request - contributions are appreciated!
Cheers, r.
participants (2)
-
hugo.levy＠gmail.com -
Robert Cimrman