I am using FEniCs to apply a load on particular point of a square mesh. When checking people has proposed to use Pointwise but it seems to be a method to fix a displacement, others has proposed to use PointSource but this also is used for one component of a vector or I don’t how to use it correctly. Any help please. Here is the code I am using: from future import print_function #!/usr/bin/env python3 -- coding: utf-8 -- “”" Created on Thu Dec 24 18:36:09 2020 @author: pierre “”" #""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""# from fenics import * from mshr import * import numpy as np Nx = 1 Ny = 1 L = 1Nx H = 1Ny mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny) E = Constant(1e9) nu = Constant(0.3) beta = E/(1-nu**2) C = beta*as_matrix([[1, nu, 0],[nu, 1, 0],[0, 0, (1-nu)/2]]) #""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""# def eps(v): return sym(grad(v)) def strain2voigt(e): “”“e is a 2nd-order tensor, returns its Voigt vectorial representation”"" return as_vector([e[0,0],e[1,1],2*e[0,1]]) def voigt2stress(s): “”" s is a stress-like vector (no 2 factor on last component) returns its tensorial representation “”" return as_tensor([[s[0],s[2]],[s[2],s[1]]]) def sigma(v): return voigt2stress(dot(C,strain2voigt(eps(v)))) #""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""# class Top(SubDomain): def inside(self, x, on_boundary): return near(x[1],L) and on_boundary class Left(SubDomain): def inside(self, x, on_boundary): return near(x[0],0) and on_boundary class Bottom(SubDomain): def inside(self, x, on_boundary): return near(x[1],0) and on_boundary class point(SubDomain): def inside(self, x, on_boundary): return near(x[0],Nx) and near(x[1],Ny) #""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""# exterior facets MeshFunction facets = MeshFunction(“size_t”, mesh, 1) facets.set_all(0) Top().mark(facets, 1) Left().mark(facets, 2) Bottom().mark(facets, 3) point().mark(facets,4) ds = Measure(‘ds’)[facets] #""" “”"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""# Define function space V = VectorFunctionSpace(mesh, ‘Lagrange’, 2) #""" “”"""""""""""""""""""""""""""""""""""""""""""""""""# Define variational problem du = TrialFunction(V) u_ = TestFunction(V) u = Function(V, name=‘Displacement’) a = inner(sigma(du), eps(u_))*dx #""" “”""""""""""""""""""""""""""""""""""""""""""""""# uniform traction on top boundary T = Constant((10,10)) l = dot(T, u_)*ds(4) #""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""# symmetry boundary conditions bc = [DirichletBC(V, Constant((0.,0.)), facets, 2), DirichletBC(V, Constant((0.,0.)), facets, 3), DirichletBC(V, T, point(), method=“pointwise”)] #""" “”"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""# solve(a == l, u, bc) #""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" import matplotlib.pyplot as plt plt.figure(1) plot(u, title = “u”) plt.show()

values[0] is of the dimension of a scalar function space (print x to see this). This means that you are evaluating your function at a set of coordinates at once. This, if you want to control a specific value in this array, you need to identify which entry you would like to change (by inspecting x).

Hi, Hope you are ok. I am trying to solve a 3D (Orthotropic) linear elasticity in FEniCs but still getting an error message that I am not understanding very well. Code from ufl import * from dolfin import * import numpy as np import matplotlib.pyplot as plt #*********************************************************************# N = 2 mesh = UnitCubeMesh(N, N, N) # plt.figure(1) # plot(mesh) # plt.show() #************************************************************# Ex, Ey, Ez, nuxy,nuxz,nuyz, Gxy, Gxz, Gyz = 100., 10.,1., 0.3, 0.2, 0.1, 5., 4., 3. S = as_matrix([[1./Ex,-nuxy/Ex,-nuxz/Ex,0.,0.,0.],[-nuxy/Ex,1./Ey,-nuyz/Ey,0.,0.,0.],[-nuxz/Ex,-nuyz/Ey,1/Ez,0.,0.,0.],[0.,0.,0.,1/Gxy,0.,0.],[0.,0.,0.,0.,1/Gxz,0.],[0.,0.,0.,0.,0.,1/Gyz]]) C = inv(S) print('S',S) print('C',C) #******************************************************************# def eps(v): return sym(grad(v)) def strain3voigt(e): """e is a 2nd-order tensor, returns its Voigt vectorial representation""" return as_vector([e[0,0],e[1,1],e[2,2],2*e[0,1],2*e[0,2],2*e[1,2]]) def voigt3stress(s): """ s is a stress-like vector (no 2 factor on last component) returns its tensorial representation """ return as_tensor([[s[0],s[3],0.],[s[3],s[1],0.],[0.,0.,s[2]]]) def sigma(v): return voigt3stress(dot(C,strain3voigt(eps(v)))) #***************************************************************************# class Top(SubDomain): def inside(self, x, on_boundary): return near(x[2],1) and on_boundary class Left(SubDomain): def inside(self, x, on_boundary): return near(x[1],0) and on_boundary class Bottom(SubDomain): def inside(self, x, on_boundary): return near(x[2],0) and on_boundary class Right(SubDomain): def inside(self, x, on_boundary): return near(x[1],1) and on_boundary class Direct(SubDomain): def inside(self, x, on_boundary): return near(x[0],0) and on_boundary #***********************************************************************# # exterior facets MeshFunction facets = MeshFunction("size_t", mesh, 2) facets.set_all(0) Top().mark(facets, 1) Left().mark(facets, 2) Bottom().mark(facets, 3) Right().mark(facets,4) Direct().mark(facets,5) ds = Measure('ds', subdomain_data=facets) #***********************************************************************# # Define function space V = VectorFunctionSpace(mesh, 'CG', 2) #**************************************************************************# # Define variational problem du = TrialFunction(V) u_ = TestFunction(V) u = Function(V, name='Displacement') a = inner(sigma(du), eps(u_))*dx #**************************************************************************# # uniform traction on top boundary T = Constant((0.,0.,1e-3)) l = dot(T, u_)*ds(1) #*********************************************************************# # symmetry boundary conditions bcs = [DirichletBC(V.sub(0), Constant(0.), facets, 2), DirichletBC(V.sub(1), Constant(0.), facets, 3), DirichletBC(V.sub(2), Constant(0.), facets, 5)] #********************************************************************# u = Function(V) problem = LinearVariationalProblem(a, l, u, bcs) solver = LinearVariationalSolver(problem) solver.solve() solver.parameters['linear_solver'] = 'cg' solver.parameters['preconditioner'] = 'ilu' cg_prm = solver.parameters['krylov_solver'] cg_prm['absolute_tolerance'] = 1E-7 cg_prm['relative_tolerance'] = 1E-4 cg_prm['maximum_iterations'] = 1000 Error message UFLException: adj_expr not implemented for dimension 6. Is it my Voigt notation that is not correct ? How to correct my program to solve a 3D orthotropic linear elasticity in FEniCs? please help me

Hi, you cannot use the UFL inv operator to invert the 6x6 compliance matrix. As said in the corresponding doc, inv is implemented only for 2x2 or 3x3 since it uses a symbolic formula for the inversion. You must therefore either implement yourself the expression of C or go through numpy for instance to invert S as a numpy array first.

Hi, I hope you are ok, I am trying to solve a linear elasticity problem. [Figure 2021-06-02 162139] But my objective is to assign different Young modulus to the different parts of the mesh above. For instance I want to assign E1 to the yellow part, E2 to the green part… Please can you show me the way to do so? I am reading documents but still not understanding very well. Thanks

See for instance: [image] How to define different materials / importet 3D geometry (gmsh) Following Dolfinx discontinous expression, this code can be modified for the syntax in dolfin as follows: import dolfin mesh = dolfin.UnitCubeMesh(9, 9, 9) V = dolfin.FunctionSpace(mesh, "DG", 0) v = dolfin.Function(V) x = V.tabulate_dof_coordinates() mf = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()) class Obstacle(dolfin.SubDomain): def inside(self, x, on_boundary): return (dolfin.between(x[1], (0.2, 0.7)) and dolfin.between(x[0], (0.3, 0.7))) Obstacl… Or https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomain

Hi Dokken, Based on your advice to see the link below: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomain I have developed the code below: from __future__ import print_function #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Jun 3 08:36:47 2021 @author: pierre """ #***********************************************************# from fenics import * from mshr import * import numpy as np #***********************************************************# # Create the mesh Nx, Ny = 1 , 1 L, H = 1, 1 mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, 'left') #************************************************************# E0 = Constant(2.1e11) nu0 = Constant(0.3) beta0 = E0/(1-(nu0*nu0)) C0 = as_matrix([[beta0, beta0*nu0, 0],[beta0*nu0, beta0, 0],[0, 0, beta0*(1-nu0)/2]]) #*************************************************************# E1 = Constant(70e9) nu1 = Constant(0.3) beta1 = E1/(1-(nu1*nu1)) C1 = as_matrix([[beta1, beta1*nu1, 0],[beta1*nu1, beta1, 0],[0, 0, beta1*(1-nu1)/2]]) #*************************************************************# def eps(v): return sym(grad(v)) def strain2voigt(e): """e is a 2nd-order tensor, returns its Voigt vectorial representation""" return as_vector([e[0,0],e[1,1],2*e[0,1]]) def voigt2stress(s): """ s is a stress-like vector (no 2 factor on last component) returns its tensorial representation """ return as_tensor([[s[0],s[2]],[s[2],s[1]]]) def sigma0(v): return voigt2stress(dot(C0,strain2voigt(eps(v)))) #*************************************************************# def sigma1(v): return voigt2stress(dot(C1,strain2voigt(eps(v)))) #*************************************************************# # Create classes for defining parts of the boundaries class Left(SubDomain): def inside(self, x, on_boundary): return near(x[0], 0.0) class Right(SubDomain): def inside(self, x, on_boundary): return near(x[0], 1.0) class Bottom(SubDomain): def inside(self, x, on_boundary): return near(x[1], 0.0) class Top(SubDomain): def inside(self, x, on_boundary): return near(x[1], 1.0) #*************************************************************# # Create classes for defining parts of the interior of the domain tol = 1E-14 class Omega_0(SubDomain): def inside(self, x, on_boundary): return x[1] = 0.5 - tol #*************************************************************# # Initialize sub-domain instances left = Left() top = Top() right = Right() bottom = Bottom() Omega_0 = Omega_0() Omega_1 = Omega_1() #*************************************************************# # Initialize mesh function for interior domains domains = MeshFunction("size_t", mesh, mesh.topology().dim()) domains.set_all(0) Omega_0.mark(domains, 1) Omega_0.mark(domains, 2) #*************************************************************# # Initialize mesh function for boundary domains boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) boundaries.set_all(0) left.mark(boundaries, 1) top.mark(boundaries, 2) right.mark(boundaries, 3) bottom.mark(boundaries, 4) #*************************************************************# # Define function space V = VectorFunctionSpace(mesh, 'CG', 1) #*************************************************************# # Define variational problem du = TrialFunction(V) u_ = TestFunction(V) u = Function(V, name = 'Displacement') #*************************************************************# # Define new measures associated with the interior domains and # exterior boundaries dx = Measure('dx', domain=mesh, subdomain_data=domains) ds = Measure('ds', domain=mesh, subdomain_data=boundaries) #*************************************************************# # Define Dirichlet boundary conditions at bottom boundary bcs = DirichletBC(V, Constant((0., 0.)), boundaries, 4) #*************************************************************# a = inner(sigma0(du), eps(u_))*dx(1) + inner(sigma1(du), eps(u_))*dx(2) #*************************************************************# # uniform traction on top boundary T = Constant((0, 1e7)) l = dot(T, u_)*ds(2) #*************************************************************# solve(a == l, u, bcs) But I have the following error message: *** Error: Unable to set given (local) rows to identity matrix. *** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal). *** Where: This error was encountered inside PETScMatrix.cpp. *** Process: 0 I tried to solve the problem but unsuccessfully. Please I need your help

Sorry Dokken, I see my error Thanks

Hi, Thank you again for the help, I have now solved the problem of defining subdomains with different Young modulus. Now suppose I have a domain meshed with 3 node triangles like the figure below: [Figure 2021-06-03 114705] And I want to give two different Young modulus to triangles among the triangles in the domain. How can I do that? I have tried with what is done in the link below: https://fenicsproject.discourse.group/t/how-to-assign-different-material-properties-to-elements-linear-elasticity/3742 But I was not able to solve my problem. Hope I make sens.

You need to make a minimal example, as I have already linked you to how to make a discontinuous dolfin function that can store different material properties on each cell. Please follow the structure of the posts highlighted in Read before posting: How do I get my question answered? - #4 by nate

How to apply a load on a single point or a node

Uncategorized

dokken June 2, 2021, 8:36pm 14

See for instance:

Or
https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomain

Topic		Replies	Views
Applying a force on a face variational formulation	25	3561	April 8, 2020
Help a starter! variational formulation	26	3524	April 8, 2020
How to assign different material properties to elements (Linear Elasticity) Uncategorized	19	3980	July 20, 2020
Point load and body force variational formulation	70	7139	March 13, 2020
How to apply many forces when solving elasticity problem dolfinx	11	829	June 4, 2025

How to apply a load on a single point or a node

Related topics