Unexpected resutls for larger system

General
1

Hello everyone,

I am running a linear elastic problem using Fenics. The mesh files are generated in Gmsh. I found that it works for smaller domain, while the results are all zeros for larger domain. It seems the displacement boundary conditions were not actually applied, but there are no errors during the running. The job finished normally, but the displacements in the output file are all zeros. Does anyone have any idea about this issue? Any comments and suggestions are appreciated.

from __future__ import print_function
from dolfin import *
from ufl import nabla_div
from mpi4py import MPI
import numpy as np
import meshio
import gmsh
import sys

if comm.rank == 0:

    gmsh.initialize()
    gmsh.merge('sample1_600.stl')

    s = gmsh.model.getEntities(2)
    l = gmsh.model.geo.addSurfaceLoop([e[1] for e in s])
    V0 = gmsh.model.geo.addVolume([l])

    gmsh.model.addPhysicalGroup(3, [V0], 1)
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(3)

    #obtain the model dimension
    dimensionBox = gmsh.model.getBoundingBox(-1, -1)
    gmsh.write("S1_600.msh")
    gmsh.finalize()

   def create_mesh(mesh, cell_type, prune_z=False):
        cells = mesh.get_cells_type(cell_type)
        cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
        points = mesh.points[:, :2] if prune_z else mesh.points
        out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
        return out_mesh

        mesh_from_file = meshio.read("S1_600.msh")
        tetra_mesh = create_mesh(mesh_from_file, "tetra", prune_z=False)
        meshio.write("tetra_mesh.xdmf", tetra_mesh)
else:
    dimensionBox = None
comm.Barrier()

dimensionBox = comm.bcast(dimensionBox, root=0)

xmin = dimensionBox[0]
ymin = dimensionBox[1]
zmin = dimensionBox[2]
xmax = dimensionBox[3]
ymax = dimensionBox[4]
zmax = dimensionBox[5]


mesh = Mesh()
with XDMFFile("tetra_mesh.xdmf") as infile:
     infile.read(mesh)

V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-6

def boundary_Left(x, on_boundary):
    return on_boundary and near(x[0], xmin, tol)

def boundary_Right(x, on_boundary):
    return on_boundary and near(x[0], xmax, tol)

def boundary_Front(x, on_boundary):
    return on_boundary and near(x[1], ymin, tol)

def boundary_Back(x, on_boundary):
    return on_boundary and near(x[1], ymax, tol)

def boundary_Top(x, on_boundary):
    return on_boundary and near(x[2], zmax, tol)

def boundary_Bottom(x, on_boundary):
    return on_boundary and near(x[2], zmin, tol)

applied_strain = 1.0e-4
applied_disp = applied_strain * (xmax - xmin)

bc_l = DirichletBC(V.sub(0), Constant(0), boundary_Left)
bc_r = DirichletBC(V.sub(0), Constant(applied_disp), boundary_Right)
bc_f = DirichletBC(V.sub(1), Constant(0), boundary_Front)
bc_ba = DirichletBC(V.sub(1), Constant(0), boundary_Back)
bc_bo = DirichletBC(V.sub(2), Constant(0), boundary_Bottom)
bc_t = DirichletBC(V.sub(2), Constant(0), boundary_Top)

bcs = [bc_l, bc_r, bc_f, bc_ba, bc_bo, bc_t]

u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))

a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx

u = Function(V)
solve(a == L, u, bcs, solver_parameters = {'linear_solver':'mumps'})

XDMFFile('disp_test1.xdmf').write(u)
2

This problem is not reproducible, and hence unlikely to receive much attention from the maintainers of the forum. Please strive to make your example reproducible.

3

Thank you for the comments. I modified the post. How should I upload the .stl input file here?

4

Add it as a plain text with 3x` encapsulation, or if it is top big, use Google drive or similar tools to share a link to it.

5

Thanks. The link to the input model file can be found as follows:
https://1drv.ms/u/s!AjzJPj4bEblbgRPn3zfHepVfIste?e=Qb1hqL

6

I figured it out what was the problem. I need to increase the tolerance for defining boundary sets. Previously, I used 1.0e-6, that was too small for larger domain. Hence, the boundary sets were not identified and the boundary conditions did not apply.

1 Like