Point sources in CBC

Hi Yaakoub,

I have not use pointsource very often, but I made this little piece of code work:

from cbc.pdesys import *
from numpy import linspace

mesh = UnitSquare(10, 10)
problem_parameters['time_integration'] = 'Steady'
problem = Problem(mesh, problem_parameters)

solver_parameters['degree']['p'] = 2
solver_parameters['degree']['c'] = 0
solver_parameters['family']['c'] = 'R'
solver_parameters['iteration_type'] = 'Picard'
Potential = PDESystem([['p', 'c']], problem, solver_parameters)

pointsources = []
for x in linspace(0.2, 0.8, 10):
    pointsources.append(PointSource(Potential.V['p'], Point(x, x), -0.01))

# Modify apply interface
PointSource.apply1 = PointSource.apply
def myapply(self, *args):
    if len(args) == 2:
        self.apply1(args[1])
    elif len(args) == 1:
        self.apply1(args[0])
    elif len(args) == 3:
        self.apply1(args[1])
PointSource.apply = myapply

class potential(PDESubSystem):
    def form(self, p, v_p, c, v_c, **kwargs):
        return inner(grad(p), grad(v_p))*dx + \
            (c*v_p + v_c*p)*dx + Constant(0)*v_p*dx

Potential.add_pdesubsystem(potential, ['p', 'c'], bcs=pointsources)
problem.solve()

I need to change the apply interface to make it work for all versions of apply(b), apply(A, b) and apply(A, b, x). Have not tested it much though, so perhaps you need to refine it.

Best regards

Mikael

Den Sep 19, 2012 kl. 7:51 PM skrev Yaakoub El Khamra:

> Question #209031 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/209031
>
> Description changed to:
> I have a general question about point/line/surface sources/sinks in CBC.
> At the moment I am doing something close to:
>
> class Qw_Expression(Expression):
> def eval(self, values, x):
> c=[0.5, 0.5]
> if (CalcDistance2D(x,c) <= DOLFIN_EPS):
> return -12.0
>
> Qw = Qw_Expression()
>
> and adding Qw*v*dx in the variational form. How does this compare
> against Point/PointSource in dolfin? Should I switch to the
> Point/PointSource class in dolfin? Do I have to do anything special to
> use it in CBC besides apply it? Also how do I apply a point source to a
> PDESys?
>
> Slight aside: I can create a similar expression to add line/surface
> sources/sinks. Is it fine if I do that or should I be doing something
> else?
>
> Regards
> Yaakoub
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.