How to compute $f<g$ amost everywhere for trialfunctions?

371
views
0
9 months ago by
I am working with finite elements for Maxwell's Equations (i.e. with Nedelec's edge elements) and for computation I'm using the <a href="https://fenicsproject.org/">FEniCS-project</a>. While implementing the Augmented Lagragian Method, I need to solve a PDE with a non-linearity of the form:

\begin{align*}
f(v) := \begin{cases}
v, \quad &\text{ on } A=\lbrace x : |v(x)| \leq g(x)\rbrace\\
\frac{v}{|v|}g, &\text{ on } \Omega \backslash A
\end{cases}
\end{align*}

I wanted to implement this by defining a function:
import dolfin

#Define initial values and data
f = dolfin.Constant((1.0,1.0,1.0))
bound = dolfin.Constant(1.0)

#Create mesh and define function spaces
nx = ny = nz = 10
mesh = dolfin.UnitCubeMesh(nx, ny, nz)
V = dolfin.FunctionSpace(mesh, 'N1curl', 1)

E = dolfin.TrialFunction(V)
v = dolfin.TestFunction(V)

#Define linear and bilinear forms
a = dolfin.dot(E, v)*dolfin.dx + dolfin.dot(dolfin.curl(E), dolfin.curl(v))*dolfin.dx
F = dolfin.dot(f, v)*dolfin.dx

#Define non-linearity
def non_lin(func, bound):
if func < bound:#somehow almost everywhere
return func
else:
return func / abs(func) * bound

a_new = a + non_lin(E, bound)

u = dolfin.Function(V)
dolfin.solve(a_new == F,u)
​

Thats an example which is not working yet... The problem lies in the function non_lin. Can someone help me?

Community: FEniCS Project
There are more problems with your code.
I don't even understand the first lines..
f = dolfin.Constant((1.0,1.0,1.0))
bound = f = dolfin.Constant(1.0)
​
what is f now?

Are you sure you want to put "E" into your non_lin function?
written 9 months ago by Lukas O.
Thanks for your answer. Actually, this was a typo. The function $f$ in my question was somehow not correct, as well. I edited it.
written 9 months ago by FFoDWindow

0
9 months ago by
Restate your problem to solve with NewtonSolver --- see the FEniCS tutorial for more info --- then use "conditional"

I've updated my function $f$ in the question. I guess with this it's not possible to use "conditional", is it? Is there any other possibility? Thanks.