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 intU
and 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.