Skip to content
Snippets Groups Projects
Commit 25580fbf authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove elasticity benchmarks of the summer school

parent f56897c6
No related branches found
No related tags found
No related merge requests found
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthFromPoints = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 1;
unit = 1.0e-02 ;
e1 = 4.5 * unit ;
e2 = 6.0 * unit / 2.0 ;
e3 = 5.0 * unit / 2.0 ;
h1 = 5.0 * unit ;
h2 = 10.0 * unit ;
h3 = 5.0 * unit ;
//h4 = 1.0 * unit ;
h4 = 2.0 * unit ;
h5 = 4.5 * unit ;
R1 = 1.0 * unit ;
R2 = 1.5 * unit ;
//R2 = 1.5 * unit ;
r = 2 * unit ;
ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
ssin = ( 1.0 - ccos*ccos )^0.5 ;
eps = 0.01 * unit;
Lc1 = 0.001 ;
Lc2 = 0.001 ;
Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};
Point(2) = { -e1-e2, h1 , 0.0 , Lc1};
Point(3) = { -e3-r , h1 , 0.0 , Lc2};
Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
Point(5) = { -e3 , h1+r , 0.0 , Lc2};
Point(6) = { -e3 , h1+h2-eps, 0.0 , Lc1};
Point(7) = { e3 , h1+h2, 0.0 , Lc1};
Point(8) = { e3 , h1+r , 0.0 , Lc2};
Point(9) = { e3+r , h1+r , 0.0 , Lc2};
Point(10)= { e3+r , h1 , 0.0 , Lc2};
Point(11)= { e1+e2, h1 , 0.0 , Lc1};
Point(12)= { e1+e2, 0.0 , 0.0 , Lc1};
Point(13)= { e2 , 0.0 , 0.0 , Lc1};
Point(14)= { R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
Point(15)= { 0.0 , h5 , 0.0 , Lc2};
Point(16)= { -R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
Point(17)= { -e2 , 0.0 , 0.0 , Lc1};
Point(18)= { -R2 , h1+h3 , 0.0 , Lc2};
Point(19)= { -R2 , h1+h3+h4, 0.0 , Lc2};
Point(20)= { 0.0 , h1+h3+h4, 0.0 , Lc2};
Point(21)= { R2 , h1+h3+h4, 0.0 , Lc2};
Point(22)= { R2 , h1+h3 , 0.0 , Lc2};
Point(23)= { 0.0 , h1+h3 , 0.0 , Lc2};
Point(24)= { 0 , h1+h3+h4+R2, 0.0 , Lc2};
Point(25)= { 0 , h1+h3-R2, 0.0 , Lc2};
Line(1) = {1 ,17}; /* ux=uy=0 */
Line(2) = {17,16};
Circle(3) = {14,15,16};
Line(4) = {14,13};
Line(5) = {13,12}; /* ux=uy=0 */
Line(6) = {12,11};
Line(7) = {11,10};
Circle(8) = { 8, 9,10};
Line(9) = { 8, 7};
Line(10) = { 7, 6}; /* T=10000 N */
Line(11) = { 6, 5};
Circle(12) = { 3, 4, 5};
Line(13) = { 3, 2};
Line(14) = { 2, 1};
Line(15) = {18,19};
Circle(16) = {21,20,24};
Circle(17) = {24,20,19};
Circle(18) = {18,23,25};
Circle(19) = {25,23,22};
Line(20) = {21,22};
Line Loop(21) = {17,-15,18,19,-20,16};
//Plane Surface(22) = {21};
Line Loop(23) = {11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10};
Plane Surface(24) = {23,21};
//Physical Line(25) = {9,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,19,20,10};
//Physical Surface(26) = {22,24};
Physical Line(9) = {11};
Physical Line(10) = {10};
Physical Line(8) = {1, 5};
Physical Surface(7) = {24};
Physical Point(20) = {2};
Physical Point(21) = {11};
// Physical Line(25) = {11, 12, 13, 14};
function neumannCondition (x,y,z)
return {-100e6,0,0}
end
function errorOnView (view, dim)
basis = polynomialBasis.find(2)
gaussPoints = fullMatrix(1,1)
gaussWeights = fullMatrix(1,1)
gaussIntegration.getTriangle(3,gaussPoints,gaussWeights)
nnodes = 3 -- only triangles
viewData = view:getData()
values = fullMatrix(nnodes,dim)
nodes = fullMatrix(nnodes,3)
gaussXYZ = fullMatrix(gaussPoints:size1(),3)
gaussValues = fullMatrix(gaussPoints:size1(),dim)
nEntities = viewData:getNumEntities(-1)
psi = fullMatrix(gaussPoints:size1(),3)
basis:f(gaussPoints,psi)
for i=0,nEntities-1 do
nElements = viewData:getNumElements(-1,i)
for j=0,nElements-1 do
if(viewData:getNumNodes(0,i,j)==nnodes) then
viewData:getAllValuesForElement(0,i,j,values)
viewData:getAllNodesForElement(0,i,j,nodes)
end
gaussXYZ:gemm(psi,nodes,1,0)
gaussValues:gemm(psi,values,1,0)
end
end
end
m = GModel()
m:load("conge.geo")
m:load("conge.msh")
e = elasticitySolver(m,1)
e:addElasticDomain (7, 200e9, 0.3)
e:addDirichletBC (1,8,0,0)
e:addDirichletBC (1,8,1,0)
e:addDirichletBC (1,8,2,0)
e:addNeumannBCLua (1,9,"neumannCondition")
sys = linearSystemCSRTaucs()
e:assemble(sys)
sys:systemSolve()
view = e:buildVonMisesView("vonMises")
view = e:buildDisplacementView("displacement")
view:write("displacement.msh",5,false)
errorOnView(view,3)
lcExt = 0.1;
lcIn = 0.01;
a=0.05;
L=3;
Point(1)={-L/2,-L/2,0,lcExt};
Point(2)={L/2,-L/2,0,lcExt};
Point(3)={L/2,L/2,0,lcExt};
Point(4)={-L/2,L/2,0,lcExt};
Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};
Point(10) = {0,0,0,lcIn};
Point(11) = {a,0,0,lcIn};
Point(12) = {0,a,0,lcIn};
Point(13) = {-a,0,0,lcIn};
Point(14) = {0,-a,0,lcIn};
Circle (10)= {11,10,12};
Circle (11)= {12,10,13};
Circle (12)= {13,10,14};
Circle (13)= {14,10,11};
Line Loop(1)={1,2,3,4};
Line Loop(2)={10,11,12,13};
Plane Surface(1) = {1,2};
Physical Line("Bottom")={1};
Physical Line("Right")={2};
Physical Line("Top")={3};
Physical Line("Left")={4};
Physical Line("Cicrcle")={10,11,12,13};
Physical Surface("Surface")={1};
Physical Point("Center") = {11};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
a = 0.05
T = 10e9
E = 200e9
nu = 0.3
-- --[[
--pull on sides
function analyticalSolution(x,y,z)
local r = math.sqrt(x*x+y*y)
local theta = math.atan2(y,x)
local ar = a/r
local ur = T*r/(2*E) * (1-nu+(1+nu)*ar^2 + math.cos(2*theta)*((1+nu)*(1-ar^4)+4*ar^2))
local utheta = - T*r/(2*E) * math.sin(2*theta)*((1+nu)*(1+ar^4)+2*(1-nu)*ar^2)
return {ur*math.cos(theta)-utheta*math.sin(theta), ur*math.sin(theta)+utheta*math.cos(theta),0}
end
-- --]]
--[[
--pull on corner
function analyticalSolution(x,y,z)
local r = math.sqrt(x*x+y*y)
local theta = math.atan2(y,x)
local ar = a/r
local ur = math.sin(2*theta)*T*r/E * (1+nu+4*ar^2-ar^4)
local utheta = math.cos(2*theta)*T*r/E * (1+nu + 2*(1-nu)*ar^2+ar^4)
return {ur*math.cos(theta)-utheta*math.sin(theta), ur*math.sin(theta)+utheta*math.cos(theta),0}
end
--]]
function analyticalSolutionX (x,y,z)
return analyticalSolution(x,y,z)[1]
end
function analyticalSolutionY (x,y,z)
return analyticalSolution(x,y,z)[2]
end
function errorOnView (view, dim, solutionFunction)
local basis = polynomialBasis.find(2)
local gaussPoints = fullMatrix(1,1)
local gaussWeights = fullMatrix(1,1)
gaussIntegration.getTriangle(3,gaussPoints,gaussWeights)
local nnodes = 3 -- only triangles
local viewData = view:getData()
local values = fullMatrix(nnodes,dim)
local nodes = fullMatrix(nnodes,3)
local gaussXYZ = fullMatrix(gaussPoints:size1(),3)
local gaussValues = fullMatrix(gaussPoints:size1(),dim)
local nEntities = viewData:getNumEntities(-1)
local psi = fullMatrix(gaussPoints:size1(),3)
local face = m:getFaceByTag(1)
local errorViewFactory = PViewFactory("error", "ElementNodeData",m, 0, 1)
basis:f(gaussPoints,psi)
local nGaussPoints = gaussPoints:size1()
local totalError = 0
local errMat = fullMatrix(3,1)
for i=0,nEntities-1 do
local nElements = viewData:getNumElements(-1,i)
for j=0,nElements-1 do
local element = face:getMeshElement(j)
if(viewData:getNumNodes(0,i,j)==nnodes) then
viewData:getAllValuesForElement(0,i,j,values)
viewData:getAllNodesForElement(0,i,j,nodes)
gaussXYZ:gemm(psi,nodes,1,0)
gaussValues:gemm(psi,values,1,0)
for k=0,nGaussPoints-1 do
local sol = solutionFunction(gaussXYZ:get(k,0),gaussXYZ:get(k,1),gaussXYZ:get(k,2))
local err = 0
for l=0,dim-1 do
err = err + (sol[l+1]-gaussValues:get(k,l))^2
end
totalError = totalError + err*gaussWeights:get(k,0)*element:getJacobianDeterminant(gaussPoints:get(k,0),gaussPoints:get(k,1),gaussPoints:get(k,2))
end
local numVertices = element:getNumVertices()
for k=0,numVertices-1 do
local vertex = element:getVertex(k)
local sol = solutionFunction(vertex:x(), vertex:y(), vertex:z())
local err = 0
for l=0,dim-1 do
err = err + (sol[l+1]-values:get(k,l))^2
end
errMat:set(k,0,err)
end
errorViewFactory:setEntry(element:getNum(),errMat)
end
end
end
errorViewFactory:createView()
return math.sqrt(totalError)
end
m = GModel()
m:load("hole.geo")
m:load("hole.msh")
e = elasticitySolver(m,1)
e:addElasticDomain (6, E, nu)
--e:addDirichletBCLua(1,5,0,"analyticalSolutionX")
--e:addDirichletBCLua(1,5,1,"analyticalSolutionY")
e:addNeumannBC (1,5,{0,0,0})
e:addDirichletBCLua(1,1,0,"analyticalSolutionX")
e:addDirichletBCLua(1,1,1,"analyticalSolutionY")
e:addDirichletBCLua(1,2,0,"analyticalSolutionX")
e:addDirichletBCLua(1,2,1,"analyticalSolutionY")
e:addDirichletBCLua(1,3,0,"analyticalSolutionX")
e:addDirichletBCLua(1,3,1,"analyticalSolutionY")
e:addDirichletBCLua(1,4,0,"analyticalSolutionX")
e:addDirichletBCLua(1,4,1,"analyticalSolutionY")
sys = linearSystemCSRTaucs()
e:assemble(sys)
sys:systemSolve()
e:postSolve()
--view = e:buildVonMisesView("vonMises")
view = e:buildDisplacementView("displacement")
print(errorOnView(view,3,analyticalSolution))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment