I’m new to Netgen/NGSolve. I’m trying to solve a few simple elasticity problems to start with. I have a few questions in this regard.
How do we define boundaries which are only a part of an edge in 2D or face in 3D (for example: not the entire bottom, top, etc. portion of the boundary)?
1a Can this be done in the Python interface based on the geometric location of the portion of the boundary surfaces where the Dirichlet BCs need to be imposed?
Is it possible to impose Dirichlet BCs at a specific node? If so, how is this done?
How do we impose Dirichlet BCs on a specific DOF of a given node?
you can split the boundary in the geometry into multiple parts, each using one segment and then defining the bc on the subpart of the edge.
2&3. Yes you can query the dirichlet dofs of a finite element space and set them to not free if they should be dirichlet dofs
freedofs = fes.FreeDofs()
freedofs.Clear(dofnr)
You can iterate over the mesh topology, as a reference see here:
and for example set all edge and vertex dofs of edges to dirichlet where both vertices have an x-coordinate smaller than 0:
for edge in mesh.edges:
if all(mesh[v].point[0] < 0 for v in edge.vertices):
for vertex in edge.vertices:
for dof in fes.GetDofNrs(vertex):
freedofs.Clear(dof)
for dof in fes.GetDofNrs(edge):
freedofs.Clear(dof)
Note, that mesh.edges only gives you NodeIDs, you have to get the node by inserting them into the mesh and then you can get the point (and then the x value of that) from a vertex node.
I’m following up on my earlier questions. I’ve attached a file with a 2D mesh (SEC_2D_lin_tri.vol) and the python script (try2.py) to set up a simple static elasticity problem in NGSolve. The solution is a 2D vector field with components, say, ux & uy. The Dirichlet BCs to be applied are:
a) uy=0 on Ysym (x>=0 and y=0)
b) ux=0 at the bottom right corner (x=10, y=0)
The top surface (Pres_load) has a non-zero Neumann condition. This is a well-posed problem and is a standard test problem in linear elasticity.
When I use the function space for the trial and test functions and set the flag dirichlet=“Ysym”, both degrees-of-freedom (ux and uy) are enforced on Ysym - this is not what I want. So I do not set the dirichlet flag while declaring the funtion space. Instead, I use lines 34 to 55 in the python script to define the Dirichlet BCs similar to your suggestion. When I print the FreeDofs() bit array, all the bits corresponding to the Dirichlet BCs are ‘0’, and I verified that the correct DOFs are selected. However, the system matrix seems to be singular and the solution is NaN everywhere.
Can you let me know what I’m getting wrong here? Is there a cleaner way to implement BCs on individual DOFs (something similar to line 92)?
Another unrelated issue: When I read this mesh and try a uniform refinement, mesh.Refine(), the code throws a segmentation fault.
The VectorH1 space already has as many components as the mesh dimension. If you set dim there you will get dim * mesh.dim dimensions, this is not what you want.
fes = VectorH1(mesh, order=1)
You can set your boundary conditions this way:
fesx, fesy = fes.components
for el in fesx.Elements(BND):
if el.mat == "Ysym":
for dof in el.dofs:
freedofs.Clear(fesx.ndof + dof) # offset of x-dofs
# find bottom right vertex
for v in mesh.vertices:
if abs(v.point[0] - 0.5*SB) + abs(v.point[1]) < tol:
for dof in fesx.GetDofNrs(v):
freedofs.Clear(dof)
Thanks for your response. I’m now able to apply Dirichlet BCs on specific DOFs at a given node.
If I try to use your method with the ‘dim’ parameter instead of VectorH1 - NGSolve throws this error " ‘components’ is available only for product spaces". What is the difference between setting VectorH1 and setting the ‘dim’ parameter? If it is too much to explain on this forum can you please point me to any NGSolve documentation if available?
Regarding mesh refinement for the previously attached .vol mesh file: as I mentioned earlier, NGSolve throws a segmentation fault when I try mesh.Refine(). Can you please let me know what the issue could be?
The difference is mainly in internal memory layout, the H1(…, dim=2) defines a space with vector valued degrees of freedom. So each dof is a vector with 2 components, each matrix entry is actually a 2x2 matrix,…
The VectorH1 is a compound space of multiple H1 spaces. Actually almost the same as doing
vech1 = FESpace([h1,h1])
with the difference, that the VectorH1 defines additional convenience functionality like taking the gradient, divergence, …
This is why it is not possible to define only the x component as dirichlet in the H1(dim) space but it is possible in the VectorH1. But one could achieve fixed boundaries for the H1 space by penalizing with robin boundary terms.
I just wanted to check back with you if you are able to duplicate the segfault issue with mesh refinement. If so, please let me know if there is a way to resolve the problem.
Hi sorry I forgot to answer.
I think I fixed that, at least in my version now I cannot reproduce a segfault. We are putting out a new release today (1907), tell me if you still segfault with it.
I have a follow-up question on this topic. Using your suggestion, I am able to apply homogeneous Dirichlet BCs on specific DOFs of a discretization of a VectorH1 space.
How can non-homogeneous Dirichlet BCs be applied on such DOFs? To be more specific, I have attached a NGSolve python file and an associated vol file. Is there a way to apply non-zero Dirichlet BCs (only in the y-direction) on the boundary marked “Pres_load” ?
First mark te appropriate dofs as Dirichlet and store them for later:
pres_load_dofs = []
for el in fesx.Elements(BND):
if el.mat == "Pres_load":
for dof in el.dofs:
freedofs.Clear(fesx.ndof+dof)
pres_load_dofs.append(fesx.ndof+dof)
Then add the non-homogeneous data to the solution grid function at the appropriate dofs and solve the system as in the documentation:
gfu.vec[:] = 0.0
for dof in pres_load_dofs:
gfu.vec[dof] = 1
# Obtain solution
res = gfu.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * res
Thank you for your detailed response. It works as expected for linear elements. But when I try this on elements of order 2, the displacements at mid-nodes are not the same as the applied value.
I print the values in gfu.vec corresponding to these DOFs and the values are as applied. But when I print the gfu(
Thank you for your detailed response. It works as expected for linear elements.
However, when I try this on elements of order 2, the displacements at mid-nodes are not the same as the applied value (uy=0.01 on the top surface). I tried to ensure that the Dirichlet BCs are applied even on the mid-nodes (edge dofs) of the top surface.
I print the values in gfu.vec corresponding to these DOFs and the values are as applied. But when I print the gfu(mesh()) values corresponding to the mid-nodes, the values are incorrect. I have attached the python script. Can you please tell me if this is just an artifact of post-processing or if there is actually an error in the solution?
Thank you for clarifying the reason for this discrepancy.
However, when I try nodal basis functions as you suggested with the following line:
fes = VectorH1(mesh=mesh, ele_order=2, type=“nodal”),
I get the following warning -
[color=blue]WARNING: kwarg ‘type’ is an undocumented flags option for class <class ‘ngsolve.comp.VectorH1’>, maybe there is a typo?, [/color]
and the error still persists on the mid-nodes of the element edges
As to your other suggestion, I am unable to use the “Set” function to apply non-homogeneous Dirichlet BCs on the hierarchical basis as there is no option in NGSolve to use the "Set " function to apply BCs on specific DOFs of vector-valued functions on the boundary. If I do, then the non-zero values are applied to all DOFs along that boundary.
When I declare the FESpace as you suggest, the code works as expected for order = 1 (this also requires other changes to the script as some utilities like div, grad, do not work for compound FE spaces - you had mentioned this earlier).
However, for any order > 1 NGSolve throws the following error:
[color=blue]catch in AssembleBilinearform 2: Inconsistent number of degrees of freedom, vb=VOL fel::GetNDof() = 12 != dnums.Size() = 6![/color]
One problem seems to be that fes.ndof is fixed at the same value (as fes.ndof for order=1) even when any order >1 is used.
I have attached the sample python script for your reference.
If there is no geometry attached to the mesh, it just takes the mid point of every edge, if the geometry is there, curves will be approximated as well.
.