Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Help - plate model

More
1 year 1 month ago #4667 by andressa
Help - plate model was created by andressa
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

 
More
1 year 1 month ago #4669 by mneunteufel
Replied by mneunteufel on topic Help - plate model
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
Time to create page: 0.165 seconds