Dear Pylith team,
Can the software impose both Dirichlet boundary conditions and Newman boundary conditions on a boundary? I want to impose both velocity boundary conditions and lithostatic pressure boundary conditions on the model boundary, but I don’t know how to implement it. Thank you in advance.
Nick
Yes, you can impose both Neumann and Dirichlet boundary conditions on the same boundary. If you impose a normal traction, then you do not want to constrain the displacement corresponding to the normal direction. Similarly, if you impose a shear traction, then you do not want to constrain the displacement corresponding to the shear direction.
Hi Brad,
Thank you for your reply.
The problem is I don’t know how to impose both Neumann and Dirichlet boundary conditions on the same boundary. When I impose Neumann boundary condition, I firstly convert the default Dirichlet BC to Neumann BC, so how to fix the dof in the same boundary. It couldn’t be better if you provide an example.
Nick
I don’t have a working full-scale example of combining Neumann and Dirichlet boundary conditions, because this situation rarely comes up in practice. However, the setup is quite simple. Here is an example of how I would create Neumann and Dirichlet boundary conditions in PyLith v3 (the same general concepts apply to previous versions of PyLith but the names and components are slightly different) on the +y boundary of a domain. Note: I have omitted the other boundary conditions to make it shorter.
bc = [bc_ypos_disp, bc_ypos_tract]
bc.bc_ypos_disp = pylith.bc.TimeDependentDirichlet
bc.bc_ypos_tract = pylith.bc.TimeDependentNeumann
# Dirichlet BC - constrain tangential DOF which is x displacement
[pylithapp.problem.bc.bc_ypos_disp]
label = boundary_ypos
constrained_dof = [0]
db_auxiliary_field = spatialdata.spatialdb.SimpleDB
db_auxiliary_field.description = Dirichlet BC +y edge
db_auxiliary_field.iohandler.filename = sheardisp.spatialdb
db_auxiliary_field.query_type = linear
# Neumann BC - normal compressive traction
[pylithapp.problem.bc.bc_ypos_tract]
label = boundary_ypos
db_auxiliary_field = spatialdata.spatialdb.UniformDB
db_auxiliary_field.description = Neumann BC +y edge
db_auxiliary_field.values = [initial_amplitude_tangential, initial_amplitude_normal]
db_auxiliary_field.data = [0.0*MPa, -5.0*MPa]
Hi Brad,
Is it possible to constrain a boundary which is not aligned to cartesian coordinate axes using Dirichlet BC in PyLith V2.2 or V3. I know that there is an up_dir variable associated with Dirichlet BC formulation in V2.2 but that does not seem to change the constraining direction.
Thank you
Viven
Current Dirichlet BC in PyLith v2 and v3 are for degrees of freedom aligned with the coordinate axes. The degrees of freedom are not tagential or normal; DOF 0=x, 1=y, 2=z.
We do have plans for v3 to allow the Dirichlet BC constraints to be provided in terms of normal and tangential to the boundary, but we still have some technical details to work out. For example, what should happen at corners of a domain that are not perpendicular.
Thank you, looking forward to that release.
Viven
The idea here is to use a transformation between the local and global bases. We have implemented this in small tests. It is just some programming to hook everything up inside of PyLith.
Thanks,
Matt
Thanks again, I hope that this will work for curved interfaces also (fig below). Let me know if you would like to run tests on this mesh.
Viven
What condition do you want on that curved interface?
My interest is in constraining the curved interface by Dirichlet BC in the surface normal direction.
Viven
I am trying to understand what this is physically modeling. For example, if you have two dissimilar materials with a curved interface, and you push on the side, It will not slide along the curved interface, but deform due to the stresses normal to the interface.
If instead you tell me that the interface is somehow very weak compared to the surrounding rock, then you should really be modeling it as a fault, rather than as a Dirichlet condition.
Right now, I cannot understand the model.
The interaction between the two materials is the research question under investigation. But I can assure you that the normal traction is significantly higher to prevent anything like fault opening with non-zero tractions type of thing from happening.
Thank you
Viven
I am not suggesting fault opening. I am asking how what you are suggesting is different from a fault, and if it is not different, why are you not modeling it as a fault?
I did try modeling as a fault but the dynamic fault model in V2.2, but that didn’t converge with the applied model parameters. The same (i.e., dynamic fault rupture model) is not yet available for V3 so I am waiting for that also. But if local Dirichlet BC can be applied, I’ll prefer using that over the dynamic fault model as the former is computationally less intensive (no non-linear solver needed).
Thank you
Viven
Following up on the nonplanar geometry issue…
Modeling only one side of a fault with Dirichlet boundary conditions can be a reasonable approximation of the problem with both sides of the fault when there is only shear motion parallel to the fault. As soon as you introduce nonplanar geometry, this will not be the case and modeling only one side of the fault will be a poor approximation of the full problem. For a meaningful solution to deformation on a fault with nonplanar geometry, you should embed the fault in the middle of your domain so that you can capture the full deformation field. The boundaries should be “far” away from the fault, where “far” is at least a few fault lengths.
Thanks for that tip, Brad. I’ll keep it in mind while setting up my model.
On a side note, could you please explain why is it that with gravity (g=9.81m/s2), the vertical stresses on the top surface of the domain (depth=0) are not equal to zero? Could you point out the governing equation from the manual for this.
In your plot of the stresses output from the simulation, the values are at the centroid of each cell. You can see how the color is uniform within each cells. As a result, there are no values exactly at the ground surface. If you plot the output values at their coordinates (for example, the cell centroids) versus the analytical solution, you will find that they match.
In PyLith v2, the finite-element basis functions are hardwired to first order (basis order = 1), so the stresses and strains are constant within triangular cells (basis order = 0) and nearly constant within quadrilateral cells (basis order = 0). For output of stress and strain on a mesh with triangular cells, there will always be 1 point per cell, which is at the centroid. For output of stress and strain on a mesh with quadrilateral cells, the default output is at the quadrature points, but most visualization tools don’t know how to handle this. Consequently, we provide CellFilterAvg, which averages the values over the quadrature points before writing the output, resulting in 1 point per cell.
In PyLith v3, the user can select the finite-element basis order for the computation as well as the finite-element basis order of the output. This means, you can use a basis order of 2 for the displacement field in a computational with gravitational body forces and set the basis order of the Cauchy stress and strain derived (output) fields to 1. In such a case, the output will be at the vertices of the cells and you will get values exactly at the ground surface.
Thank you for the explanation.
Viven
