Point sources in CBC

Asked by Yaakoub El Khamra on 2012-09-19

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

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:

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.

Yaakoub El Khamra (yelkhamra) said : #2

Mikael: Thank you very much sir!

Yaakoub El Khamra (yelkhamra) said : #3

Thanks Mikael Mortensen, that solved my question.

Yaakoub El Khamra (yelkhamra) said : #4

So I just tried this code again with the latest dolfin from development and it breaks with the following error:

---
*** Error: Unable to create intersection operator.
*** Reason: IntersectionOperatorImplementation is not available, DOLFIN has been compiled without CGAL.
*** Where: This error was encountered inside IntersectionOperator.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

Looking at the cmake for dolfin it was compiled with CGAL. Can you please confirm that this also breaks for you and that it is indeed a dolfin issue?

Many thanks in advance.

Regards
Yaakoub