dolfin team mailing list archive

[Neilen Marais <question152565@xxxxxxxxxxxxxxxxxxxxx>]

Question #152565 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/152565

    Status: Answered => Solved

Neilen Marais confirmed that the question is solved:
I implemented my own point-source that works with Nédélec elements, see
code below. I'm not 100% sure if it is mathematically correct for the
case where the point-source is co-incident with a cell boundary (since
I'm only considering one element), but seems to work OK. I'm also not
sure if it will work for other element types, but presumably you could
use dolfin.PointSource that is now included in dolfin 0.9.11. For a full
example see http://blogs.sun.ac.za/cemagg/2011/05/16/electrical-dipole-
and-first-order-absorbing-boundary-condition-abc-with-dolfin/


def calc_pointsource_contrib(V, source_coords, source_value):
    """Calculate the RHS contribution of a current point source (i.e. electric dipole)
    Input Values
    -------------
    @param V: dolfin FunctionSpace object
    @param source_coords: length 3 array with x,y,z coordinates of point source
    @param source_value: length 3 array with x,y,z componets of source current 
    Return Values
    -------------
    (dofnos, rhs_contribs) with

    dofnos -- Array of degree of freedom indices of the source
contribution

    rhs_contribs -- Numerical values of RHS contribution, such that
        RHS[dofnos] += rhs_contribs will add the current source to the system.

    """
    source_coords = N.asarray(source_coords, dtype=N.float64)
    source_value = N.asarray(source_value)
    dm = V.dofmap()
    dofnos = N.zeros(dm.max_cell_dimension(), dtype=N.uint32)
    source_pt = dolfin.Point(*source_coords)
    try:
        cell_index = V.mesh().any_intersected_entity(source_pt)
    except StandardError:
        # CGAL as used by dolfin to implement intersection searches
        # seems to break with 1-element meshes
        if dolfin.Cell(V.mesh(), 0).intersects(source_pt):
            cell_index = 0
        else: raise
    c = dolfin.Cell(V.mesh(), cell_index)
    # Check that the source point is in this element    
    assert(c.intersects_exactly(source_pt)) 

    dm.tabulate_dofs(dofnos,  c)
    finite_element = V.dolfin_element()
    no_basis_fns = finite_element.space_dimension()
    # Vector valued elements have rank of 1
    assert(finite_element.value_rank() == 1)
    # Vector valued elements have only one rank (i.e. 0) along which
    # dimensions are defined. This is the dimension that the basis
    # function value vector is. Since we have 3D Nedelec elements here
    # this should be 3
    bf_value_dimension = finite_element.value_dimension(0)
    el_basis_vals = N.zeros((no_basis_fns, bf_value_dimension), dtype=N.float64)
    finite_element.evaluate_basis_all(el_basis_vals, source_coords, c)
    # Compute dot product of basis function values and source value
    rhs_contribs = N.sum(el_basis_vals*source_value, axis=1)
    return dofnos, rhs_contribs

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.