-
Notifications
You must be signed in to change notification settings - Fork 94
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Provide a general concept of orientation that will be general for whatever dimension or polytope #51
Comments
@santiagobadia, continuing the discussion about the normals. my goal is to construct I want to achieve an API like this one: using Gridap
import Gridap: ∇
ufun(x) = VectorValue(x[1]*x[2], x[2])
grad_ufun(x) = TensorValue(x[2],x[1],0.0,1.0)
∇(::typeof(ufun)) = grad_ufun
model = CartesianDiscreteModel(partition=(4,4,4))
trian = Triangulation(model)
trian_Γ = BoundaryTriangulation(model, "boundary") # An integration mesh of the entire boundary
quad = CellQuadrature(trian, order=2)
quad_Γ = CellQuadrature(trian_Γ, order=2)
u = CellField(trian, ufun)
u_Γ = restrict(u, trian_Γ)
n = NormalVector(trian_Γ) # This is the major part missing
# Check that divergence theorem holds
@assert sum( integrate( div(u) , trian, quad) ) ≈ sum( integrate( u_Γ * n , trian_Γ, quad_Γ) )
The function My first guess is to construct the normal from the Jacobian of the CellGeomap of the boundary triangulation. I.e.: phi = CellGeomap(trian_Γ)
jac = ∇(phi) The Do you understand what I mean? Do you have perhaps a better way to compute the normals? |
Hi @fverdugo , I think that the general way to define outward normals is: 1) Define the outward normal in the reference FE, using the nullspace of the squared matrix of tangent vectors, etc; 2) Since the tangent vectors transform with the local change of basis with the Jacobian matrix What do you think? |
Yes, we can follow this approach. It is simple, and we already have all ingredients. I have just a question. This approach assumes that each face in the surface mesh touches a volume cell in some volumetric mesh (this is actually the case of our |
@fverdugo, I think that if you talk about general manifolds in arbitrary dimension, the situation is more complicated. First, you need to define the concept of orientable (or oriented, not sure) manifold. Not all manifolds are orientable. Assuming that you have an orientable manifold, you could probably define the normal component of a facet in the reference domain (in dimension+1, normal = (0,...0,1) and use the tranpose of the Jacobian to determine the normal component of the manifold. Again, given a oriented manifold, the normal would always point to the same side if some restrictions apply. It seems that one way to define an orientable map is via triangulation, enforcing all triangles be clockwise (or counterclockwise). If we assume that the manifold is the boundary of domain, the manifold is orientable and we can make the normal to point outwards using what I proposed. By the way, STL meshes should provide a consistent normal, so in this case, I think that there is no problem, it is just data. They are (as far as I can say) a representation of an orientable manifold. Not sure about STL for Mobius strips. Surely life is simpler in 3D using the vector product (and it can also be used in 2D by appending the Z-component). But there is no extension of the vector product in higher dimension. We can probably consider a particular case for 2D/3D manifolds in 3D/4D when needed. Or study a little bit more of differential geometry and try to provide an abstract definition. My differential geometry knowledge is quite basic. |
@santiagobadia work done! using Gridap
tags = [23,24,25]
model = CartesianDiscreteModel(partition=(3,4,2))
btrian = BoundaryTriangulation(model,tags)
n = NormalVector(btrian) # New function!
bquad = CellQuadrature(btrian,order=2)
q = coordinates(bquad)
n_q = evaluate(n,q)
bphi = CellGeomap(btrian)
x = evaluate(bphi,q)
writevtk(btrian,"btrian")
writevtk(x,"x",pointdata=["n"=>n_q]) |
Closed via PR #60 |
Place to discuss on this topic
The text was updated successfully, but these errors were encountered: