I want to solve the following weak form with well-defined BC:
\int_\Omega \nabla u \, A \cdot \nabla v \,d\Omega = \int_\Omega v^T f \,d\Omega, \quad with \, \Omega \in \mathbb{R}^n
where A is an elliptic matrix of dimension n\times n. The resulting FEM equation for each element normally reads:
K_e u_e = B_e^T A B_e u_e = f_e
where B_e defines derivatives of the shape functions. As far as I know, this is easy to assembly with fenics. But my new problem needs to assembly the following:
B_e^T A_e B_e u_e = f_e + B_e^T g_e
where A_e, f_e, g_e are given for each element. Is it possible to this in fenics?
Hello, what is exactly g_e ? From a mechanical point of view it looks like a pre-stress. If so, then just define it as an expression sig0 = Expression(...) or a Constant if it is the case and just include it in the variational inear form:
It is not a pre-stress. ge is a value/vector I want to add to the right hand side of the equation system. The problem is, that my values are defined on the element level and not on a global level. So I have 2 problems I cant solve.
How to add BTege for each element to fe and how to resolve the problem with different Ae for each element.
If it is element wise, it could be represented as a DG0 function?
I would assume something similar to this could work Dolfinx discontinous expression (note that this is dolfinx syntax, Which differs slightly from dolfin syntax)
Oh I am very sorry. I guess I did a false description of my problem. Let me correct my self. For instance I have the following code:
import numpy as np
from dolfin import *
L = 2;
#V = VectorFunctionSpace(mesh, 'Lagrange', 1)
mesh = IntervalMesh(1, 0, L)
VE = VectorElement("CG", interval, 1, dim=2)
V = FunctionSpace(mesh, VE)
mesh = Mesh(mesh)
# Define Dirichlet boundary
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)
# Solution variables
dx = Measure("dx")
du = TrialFunction(V)
U = Function(V)
u = U.vector()
v = TestFunction(V)
f = Constant((0, 0))
T = Constant((1, 1))
A = 1
a = inner(A*grad(du), grad(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
# Assembly system
K, f_vec = assemble_system(a, L, bc)
# Solve system
solve(K, u, f_vec)
In this case I want to change two different things. My number of nodes is 2 with DoF 2 in each node. Therefore I expect the assembly method to loop over the two nodes (in more Dimensions it would be over the number of GP I guess).
During the assembly over the nodes I want to change the system such that in each node he is using a different A and adding a different B_e^T*g_e to the right side.
As long as A is an element-wise constant, i would seriously recommend you to represent it as a DG0 function. This is the far-most simplest and most generic way of doing it. Having an VectorFunctionSpace(mesh,"DG",0) is a trivial extension.
Just to clearify and maybe to get some other usefull advices. My goal is to implement a data driven solver introduced by Ortiz. The algorithm of the solver looks like this:
.
