Fail to use assemble_system


83
views
0
22 days ago by
Hi everyone,

I want to solve the equation 
$\text{ ∇·v=δ(x) in Ω  }$ ∇·v=δ(x) in Ω    
v= −∇u in Ω 
u=0 on ∂Ω 

Hence, I need to apply PointSource function and before that I need to use assemble_system.

Herewith the part of my code

from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import numpy as np


center=np.array([0,0])

#need to make sure (0,0) is not the mesh point
x0=50
y0=25
nx=100
ny=50

mesh=RectangleMesh(Point(-x0,-y0),Point(x0,y0),nx,ny)

u_D = Constant(0)

def boundary(x, on_boundary):
return on_boundary

# Define the vector space
P1=VectorElement('P',triangle,1)
P2=FiniteElement('P',triangle,1)
P=MixedElement([P1,P2])
W=FunctionSpace(mesh,P)

u=TrialFunction(W)
u_1,u_2=split(u)
(v_1,v_2)=TestFunctions(W)

F=inner(u_1,v_1)*dx-u_2*div(v_1)*dx+div(u_1)*v_2*dx
a,L=lhs(F),rhs(F)
bc=DirichletBC(W.sub(1), 0, boundary)
A,b=assemble_system(a,L,bc)
delta=PointSource(W,Point(0,0),1)
delta.apply(b)

u=Function(W)
solve(A,u.vector(),b) ​

Then things went wrong with assemble_system function.
The error i get is 

Traceback (most recent call last):
File "/usr/lib/python3.5/code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 361, in assemble_system
b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 67, in _create_dolfin_form
function_spaces=function_spaces)
File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 77, in __init__
"(got {}).".format(len(sd)))
File "/usr/lib/python3/dist-packages/dolfin/cpp/common.py", line 2808, in dolfin_error
return _common.dolfin_error(location, task, reason)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to creating dolfin.Form.
*** Reason: Expecting subdomain data for exactly one domain(got 0)..
*** Where: This error was encountered inside form.py.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Thanks in advance for helping and best wishes!
Alice
Community: FEniCS Project

1 Answer


0
21 days ago by
Hello,

I would suggest the following changes on your code, the usage of the assemble_system with the derivative and the du as a trialfunction ... 

from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import numpy as np


center=np.array([0,0])

#need to make sure (0,0) is not the mesh point
x0=50
y0=25
nx=100
ny=50

mesh=RectangleMesh(Point(-x0,-y0),Point(x0,y0),nx,ny)

u_D = Constant(0)

def boundary(x, on_boundary):
    return on_boundary

# Define the vector space
P1=VectorElement('P',triangle,1)
P2=FiniteElement('P',triangle,1)
P=MixedElement([P1,P2])
W=FunctionSpace(mesh,P)

u= Function(W)
du=TrialFunction(W)

u_1,u_2=split(u)
(v_1,v_2)=TestFunctions(W)

F=inner(u_1,v_1)*dx-u_2*div(v_1)*dx+div(u_1)*v_2*dx

J = derivative(F, u, du)

bc=DirichletBC(W.sub(1), 0, boundary)
A,b=assemble_system(J, F,bc)
delta=PointSource(W,Point(0,0),1)
delta.apply(b)

u=Function(W)
solve(A,u.vector(),b)
​

Hope it helps,

Regards, Leonardo
Thanks a lot but why u is changed to Function(W) and du to be TrialFunction?
written 21 days ago by Alice Peng  
what i see is that your rhs(F) is returning an empty form because all of your weak formulation depends on the u function itself, that might be the reason for your problem when assembling the system from your original question...

Taking also into account that your lhs(F) is nonlinear, you need to provide the Jacobian in order to solve using Newton's method. For calculating the jacobian automatically using the derivative method, the displacement function du as a trialfunction is needed. I guess when defining the weak formulation, it is ok to use as a Function when providing the trial function via the jacobian...

In the hyperelasticity demo this is explained in details: https://fenicsproject.org/olddocs/dolfin/2016.1.0/cpp/demo/documented/hyperelasticity/cpp/documentation.html
written 21 days ago by lhdamiani  
okay, thanks but I have found an example which is quite similar to my situation, and it did work with my problem even with assemble_system. Here with the link:
https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/mixed-poisson-dual/python/documentation.html
written 20 days ago by Alice Peng  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »