Split a 2d solution to add different constants to each component and write to file

General
1

Hi there, I am trying to solve a two body problem and as a first step started with two parallel membranes.

I define the space

V=FunctionSpace(mesh, VectorElement(“P”, triangle, degree = 1, dim = 2))

Define forces on each membrane

f1 = Constant(-40.0)
f2 = Constant(40.0)

Define boundary conditions

u_D = Constant((0.,0.))

def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, u_D, boundary)

Split in two components the function and the test function (this I will need when putting unilateral conditions)

u = Function(V)
u1,u2 = split(u)
v = TestFunction(V)
v1,v2 = split(v)

Provide my variational formulation

F = dot((grad(u1)), (grad(v1)))dx+dot((grad(u2)), (grad(v2)))dx- (f1v1+f2v2)*dx

And solve

solve(F == 0, u, bc)

Now I wanto to write to one file u1+L and to one another u2-L and all my attempts failed. I can show on screen by using

plot(u1+L, title=‘Deflection1’)
plt.show()
plot(u2-L, title=‘Deflection2’)
plt.show()

To be honest, I get a warning that “Object cannot be plotted directly, projecting to piecewise linears”.
But further, all my attempts to write to file failed. I try

file1 = File(“place/u1.pvd”)
file2 = File(“place/u2.pvd”)
file1<< u1+L
file2<< u2-L

But I get "AttributeError: ‘Sum’ object has no attribute ‘_cpp_object’ "

I tried to interpolate and project the constant L before adding to u1 or u2, but I am doing something wrong and do not know what to try next.

This is probably very basic stuff, sorry for bothering you with this issue, but any help will be appreciated.

2

Please make a minimal reproducible example, as L is not defined in your snippets above, and your error is thus kot reproducible.
Also make sure you use 3x` encapsulation, ie

```python
#add code here
```
3

I moved the lines that dont work in the end


from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# Create mesh and define function space
mesh = UnitSquareMesh(20, 20)
V=FunctionSpace(mesh, VectorElement("P", triangle, degree = 1, dim = 2))


# Define loads

f1 = Constant(-40.0)
f2 = Constant(40.0)


# Define boundary condition
u_D = Constant((0.,0.))

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)


# Define variational problem
u = Function(V)  # Note: not TrialFunction!
u1,u2 = split(u)
v = TestFunction(V)
v1,v2 = split(v)


mu1=1.0
mu2=1.0

F = mu1*dot((grad(u1)), (grad(v1)))*dx+mu2*dot((grad(u2)), (grad(v2)))*dx- (f1*v1+f2*v2)*dx


# Compute solution
solve(F == 0, u, bc)

(u1,u2)=u.split()

ele=1.

plot(u1+ele, title='Deflection1')
plt.show()
plot(u2-ele, title='Deflection2')
plt.show()

file1 = File("sofo/u1.pvd")
file2 = File("sofo/u2.pvd")


#file1<< u1+ele
#file2<< u2-ele

# Curve plot along x = 0.5 comparing gap between membranes
tol = 0.001  # avoid hitting points outside the domain
y = np.linspace(0 + tol, 1 - tol, 101)
points = [(0.5, y_) for y_ in y]  # 2D points
teta_line = np.array([2*ele+u1(point)-u2(point) for point in points])
#p_line = np.array([p(point) for point in points])
plt.plot(y, teta_line, 'k', linewidth=2)
#plt.plot(y, p_line, 'b--', linewidth=2)
plt.grid(True)
plt.xlabel('$y$')
plt.legend(['gap'], loc='upper left')
#plt.savefig('poisson_nonlinear/curves.pdf')
plt.savefig('sofo/curves.png')
plt.show()




file1<< u1+ele
file2<< u2-ele
4

As ele is a floating point vale, you can simply call
u1.vector()[:] += 1.
u2.vector()[:] -= 1.
and then write the corresponding function u1 or u2 to file.

1 Like
5

Thank you. It kinda works. But I am still doing something wrong. For example how can it be that if I modify just one, only one of the two components, let’s say u2, then u1 is automatically modified as well. This should not happen, right?

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# Create mesh and define function space
mesh = UnitSquareMesh(20, 20)
V=FunctionSpace(mesh, VectorElement("P", triangle, degree = 1, dim = 2))


# Define loads

f1 = Constant(-40.0)
f2 = Constant(40.0)


# Define boundary condition
u_D = Constant((0.,0.))

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)


# Define variational problem
u = Function(V)  # Note: not TrialFunction!
u1,u2 = split(u)
v = TestFunction(V)
v1,v2 = split(v)


mu1=1.0
mu2=1.0

F = mu1*dot((grad(u1)), (grad(v1)))*dx+mu2*dot((grad(u2)), (grad(v2)))*dx- (f1*v1+f2*v2)*dx


# Compute solution
solve(F == 0, u, bc)

(u1,u2)=u.split()




#u1,u2=u.split()

print("u1 value at (0,0): %.2f" %(u1(0.,0.)))
print("u2 value at (0,0): %.2f" %(u2(0.,0.)))
#print(np.array([u1(0.,0.)]))
u2.vector()[:]+=1.  #I JUST MODIFY THIS ONE

print("u1 value at (0,0): %.2f" %(u1(0.,0.)))#THIS ONE GETS MODIFIED AS WELL
print("u2 value at (0,0): %.2f" %(u2(0.,0.)))




ele=1.

plot(u1, title='Deflection1')
plt.show()
plot(u2, title='Deflection2')
plt.show()

file1 = File("sofo/u1.pvd")
file2 = File("sofo/u2.pvd")


file1<< u1
file2<< u2

# Curve plot along x = 0 comparing p and w
tol = 0.001  # avoid hitting points outside the domain
y = np.linspace(0 + tol, 1 - tol, 101)
points = [(0.5, y_) for y_ in y]  # 2D points
u1_line = np.array([u1(point) for point in points])
u2_line = np.array([u2(point) for point in points])
plt.plot(y, u1_line, 'k', linewidth=2)
plt.plot(y, u2_line, 'b--', linewidth=2)
plt.grid(True)
plt.xlabel('$y$')
plt.legend(['u1 line','u2 line'], loc='upper left')
plt.savefig('sofo/curves.png')
plt.show()

Thank you again

6

Use deepcopy=True in the split, as explained in: Bitbucket

1 Like
7

WOW, thank you! It’s perfect now.