Help - plate model

Dear all,

I am working with the Reissner-minlin plate model, the trial source term f is given by the two-dimensional Dirac delta
function. However, the program returns an error using CoefficientFunction: 

f = LinearForm(source *du* dx).Assemble()
netgen.libngpy._meshing.NgException: Linearform must not have TrialFunction

Code:
##########################################################################################################
from ngsolve import *
from ngsolve.meshes import MakeStructured2DMesh
cx, cy = [], []
# one ball
cx.append(0.66666667)
cy.append(0.66666667)
# two ball
cx.append(0.5)
cy.append(0.25)
def gdirac(cx, cy):
    # eps = 7.499999999999972e-01 #0.1*meshsize
    eps = 2.499999999999994e-01
    coef = 1.0/(2.0*pi*eps*eps)
    dif = -((x-cx)*(x-cx) + (y-cy)*(y-cy) )/(2*eps)
    gdirX = -coef * exp(dif) * (x-cx)
    gdirY = -coef * exp(dif) * (y-cy)
    return gdirX, gdirY
E, nu, k = 10.92, 0.3, 5/6
mu  = E/(2*(1+nu))
lam = E*nu/(1-nu**2)
#shearing modulus
G = E/(2*(1+nu))
#thickness, shear locking with t=0.1
t = 0.1
R = 5
#force
order = 1
Db = E*t**3/(12*(1-nu**2))
mesh = MakeStructured2DMesh(quads=False, nx=6, ny=6)
moment=[]
# 1: cavitation
coef = (2.0*mu + 3*lam)
beta = 1.0e-10
m11, m12 = beta*coef, 0
m21, m22 = 0, beta*coef
# moment tensor: cisalhamento
moment.append(m11, m12], [m21, m22)

clamped = True
if clamped:
    fesB = VectorH1(mesh, order=order, dirichletx="circ|left", dirichlety="circ|bottom", autoupdate=True)
else:
    fesB = VectorH1(mesh, order=order, dirichletx="left", dirichlety="bottom", autoupdate=True)

fesW = H1(mesh, order=order, dirichlet="circ", autoupdate=True)
fes = FESpace( [fesB, fesW], autoupdate=True )
(beta,u), (dbeta,du) = fes.TnT()

fes.ndof  
nz = fes.ndof
Z = GridFunction(fes)
Gamma = HCurl(mesh,order=order, autoupdate=True)
with TaskManager():
    a = BilinearForm(fes)
    a += t**3/12*(2*mu*InnerProduct(Sym(grad(beta)),Sym(grad(dbeta))) + lam*div(beta)*div(dbeta))*dx
    a += t*k*G*InnerProduct(Interpolate(grad(u)-beta,Gamma), Interpolate(grad(du)-dbeta,Gamma))*dx
    a.Assemble()
print(a)
cx_, cy_ = cx[0], cy[0]
gdirX, gdirY = gdirac(cx_, cy_)
sorX = (beta*coef)*gdirX
sorY = (beta*coef)*gdirY
source = CoefficientFunction((sorX, sorY))
f = LinearForm(source *du* dx).Assemble()
# res = f.vec.CreateVector()
# res.data += f.vec - a.mat * Z.vec
# Z.vec.data += a.mat.Inverse(freedofs=fes.FreeDofs()) * res

##################################################################################################
Regards,
Andressa

 

Hi Andressa,

in sorX and sorY you use beta, which is a TrialFunction.
In Linearforms, however, only TestFunctions are allowed, therefore the error message.
If you use a CoefficientFunction not including some Trialfunctions this error should be gone.

The final expression in the LinearForm has to be a scalar. Currently source is a vector and du a scalar which will give you an error. If you make source scalar valued your file should work.

Best,
Michael