Tolerances in load boundary condition definitions

[itk22]

Hello,

I'm currently working on implementing point loads in FEM simulations using your library and have some questions about the tolerance (atol) used for defining load boundary conditions. I noticed that in some examples, including applications/fem/top_opt/box.py, the load location function uses certain constants multiplied by the domain dimensions:

Lx, Ly, Lz = 2., 0.5, 1.
Nx, Ny, Nz = 80, 20, 40

meshio_mesh = box_mesh(Nx, Ny, Nz, Lx, Ly, Lz, data_path)
jax_mesh = Mesh(meshio_mesh.points, meshio_mesh.cells_dict['hexahedron'])

def fixed_location(point):
    return np.isclose(point[0], 0., atol=1e-5)

# Note the tolerance here:    
def load_location(point):
    return np.logical_and(np.isclose(point[0], Lx, atol=1e-5), np.isclose(point[2], 0., atol=0.1*Lz+1e-5))

def dirichlet_val(point):
    return 0.

def neumann_val(point):
    return np.array([0., 0., -1e6])

From my understanding, there are two factors to consider when determining this tolerance. First, we need the tolerance (atol = c*Lz) to be larger than the mesh edge length (dz = Lz/Nz) to ensure that at least one node falls within this range. However, I also understand that the multiplication factor c should ideally be decoupled from the mesh resolution to ensure consistent results across different mesh sizes.

Could you provide some guidance on choosing suitable tolerances and discuss any recommended practices or considerations? For instance, is there an optimal range or formula to calculate c in relation to the mesh size or problem scale?

Also, I'm curious if there might be a more systematic way to handle this tolerance selection within the library itself, such as by providing a helper function to calculate suitable tolerances based on mesh size and problem scale. I believe this would greatly assist users in correctly applying loads and constraints to their simulations.

Thank you in advance for your help and guidance. I appreciate your work on this library and look forward to your response.

[tianjuxue]

Thanks for your interest. We purposely let users handle atol (rather than define it in a systematic way) since this is really problem dependent, and different users may have their own needs.
For this topology optimization problem, we use 10% (0.1*Lz+1e-5) of the boundary length to approximate a point load. Note that Neumann boundary condition (or traction force B.C.) should be applied on at least one face of an element, not on a single point. Sometimes, we need to adjust that to 5%, for example, so it is better to leave it there and let users figure out what they need.
A good practice may be to return the number of nodes (for Dirichlet condition) or number of faces (for Neumann condition) and print this information out to inform users. We will add this feature in our future version.

[itk22]

Thank you for the quick response. You mentioned that setting atol to 10% approximates a concentrated force/ point load. What is the closest approximation one could get with this approach? Would it correspond to a single element face with atol spanning 2 or 3 neighboring mesh nodes depending on the chosen element type? I agree that it might be a good idea to return the number of nodes/ faces for boundary condition, or even assert that the number is nonzero.

[tianjuxue]

I have just implemented a function to print out boundary info. Please use this to check if desired boundary conditions are imposed successfully. In you application code, problem.print_BC_info() will print some useful information.

jax-am/jax_am/fem/core.py

Lines 689 to 716 in 1d34e86

	def print_BC_info(self):
	"""Print boundary condition information for debugging purposes.
	"""
	if hasattr(self, 'neumann_boundary_inds_list'):
	print(f"\n\n### Neumann B.C. is specified")
	for i in range(len(self.neumann_boundary_inds_list)):
	print(f"\nNeumann Boundary part {i + 1} information:")
	print(self.neumann_boundary_inds_list[i])
	print(f"Array.shape = (num_selected_faces, 2) = {self.neumann_boundary_inds_list[i].shape}")
	print(f"Interpretation:")
	print(f" Array[i, 0] returns the global cell index of the ith selected face")
	print(f" Array[i, 1] returns the local face index of the ith selected face")
	else:
	print(f"\n\n### No Neumann B.C. found.")
	if len(self.node_inds_list) != 0:
	print(f"\n\n### Dirichlet B.C. is specified")
	for i in range(len(self.node_inds_list)):
	print(f"\nDirichlet Boundary part {i + 1} information:")
	bc_array = onp.stack([self.node_inds_list[i], self.vec_inds_list[i], self.vals_list[i]]).T
	print(bc_array)
	print(f"Array.shape = (num_selected_dofs, 3) = {bc_array.shape}")
	print(f"Interpretation:")
	print(f" Array[i, 0] returns the node index of the ith selected dof")
	print(f" Array[i, 1] returns the vec index of the ith selected dof")
	print(f" Array[i, 2] returns the value assigned to ith selected dof")
	else:
	print(f"\n\n### No Dirichlet B.C. found.")

As for how much atol is suitable to approximate a point load. I guess there is really no universal answer. You could do just one element face, which is perfectly legal.