Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.
Calculation of normal vectors on surface
Posted 22.12.2010, 10:10 GMT-5 Geometry, Materials, Mesh, Parameters, Variables, & Functions 7 Replies
Please login with a confirmed email address before reporting spam
I need to calculate normal vectors on a 3D geometry surface to calculate boundary conditions from that (it's a rather complex matlab code). I took the postinterp help to get a starting point and tried to calculate the normal vectors in the center of each element. I used postinterp and tried to determine nx, ny and nz, but, whatever I do, I only get NaNs. In the postinterp help, a parameter s was used to do that, this also works fine in 3D, but I am not sure how I would transform certain points to that parameter... Here's my code:
% Create sphere surface geometry
clear fem
fem.geom=geomcoerce('face',sphere3(1,'pos',[0,0,0]));
geomplot(fem,'edgelabels','on'), axis equal
% Generate mesh
fem.mesh = meshinit(fem);
fem.shape = 2;
fem.xmesh = meshextend(fem);
% Get element triangulation
a=get(fem.mesh,'el');
tri=a{3}.elem;
p=fem.mesh.p;
% Calculate center point of each triangle
pn=(p(:,tri(1,:))+p(:,tri(2,:))+p(:,tri(3,:)))/3;
% Calculate normal vectors
nx=postinterp(fem,'nx',pn)
% this gives only NaNs, why?
Any ideas? Also the interpolation in the mesh points does not work,
nx=postinterp(fem,'nx',fem.mesh.p)
gives me the same NaNs...
Nils
Please login with a confirmed email address before reporting spam
How do ou define a normal vector on the elements ?
in COMSOL you have the normal from the boundaries (nx ny nz) either in the material or the spatial frame etc, check the doc
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
I know I have the nx, ny, nz but I just can't get the calculation of it (postinterp, poseval) on an arbitary surface point working. The doc doesn't help, it just gives examples based on a parameterization which I will not be able to use for my final gemeotry.
Please login with a confirmed email address before reporting spam
but (nx,ny,nz) are defined per element on a boundary (pointing outward w.r.t. the domain) so why not define an average operator on your boundary and define Nx = aveop1(nx) ... then you would have something, in average, OK if your boundary is rather smooth and of limited extend
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
As you said - (nx,ny,nz) are defined per element. But how can I get this information per element?
Please login with a confirmed email address before reporting spam
try an arrow plot on boundaries with (nx,ny,nz) then in matlab you can retrieve these variables as any others
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
COMSOL uses "topological" or "surface" coordinates rather than actual coordinates to calculate the normal vectors. In 2D the topological coordinate is one-dimensional (since it is an edge) and in 3D it is 2-dimensional (since it is a surface).
In my opinion, the best way to calculate normal vectors on the surface seems to be a combination of the following commands:
1. flgeomud - Get up-down sub-domain numbering of geometry faces
2. flgeomfd - Get coordinates and derivatives for geometry faces
3. flgeomfn - Get normals from face derivatives
!! Important !!
POSTINTERP command works well for obtaining normal vector information on edges, however returns "NaN" for many face-normal vectors in 3D simulations. Hence I do not recommend POSTINTERP.
Now the philosophy behind calculating normal vectors :
Let number of subdomains be nd. The empty space outside the model is numbered 0.
Commands "flgeomfd" and "flgeomfn" can easily provide us with the face-normal information. But when we talk of normal vector on a surface it becomes mandatory to know how it is directed. For example if domains 1 and 2 share a face number 3, and if the normal vector say [1, 0, 0] on this surface is calculated from the the commands "flgeomfd" and "flgeomfn", we still do not have the information as to whether this normal is directed from domain 1 to domain 2 or domain 2 to domain 1. In other words, the normal vector information provided by these 2 commands lacks 'direction' information.
This is where the command "flgeomud" comes handy. flgeomud is a matrix of size (nd+1) x nf where nf = number of faces. Thus the information from a geometry g obtained using flgeomud looks like this:
>updown = flgeomud (g)
1 1 2 ... ( This row is for "up" domain )
0 0 1 ... ( This row is for "down" domain )
This means :
a) the normal of the face 1 is directed from sub-domain 0 (empty space) to the sub-domain 1
b) the normal of the face 2 is directed from sub-domain 0 (empty space) to the sub-domain 1
c) the normal of the face 3 is directed from sub-domain 1 to the sub-domain 2
Thus the normal vector [1,0,0] of the shared face 3 we obtained earlier points out "from" domain 1 "to" domain 2. Thus it is an "outward" pointing normal for face 1 when domain 1 is concerned and "inward" pointing normal when domain 2 is concerned.
To illustrate this I shall give an example here which presumes that there is only only one subdomain. Thus, we have
empty space : 0
subdomain : 1 = nd
A simple example is a parallelepiped (cuboid) which has six faces. This geometry has faces that are "shared" with the empty space only. Let us think of a situation where there is some heat generated inside this box and that we want to calculate the total heat energy that escapes through all the walls. For this we need to find all the face normal vectors which point "outward" from the box "into" the empty space.
==============================
% Create a simple cubic block
>b=block3(10,20,30);
>clear s
s.objs={b};
s.name={'Block'};
s.tags={'Block'};
fem.draw=struct('s',s);
fem.geom=geomcsg(fem);
% Obtain number of faces (boundaries)
>nb = (geominfo(b, 'out',{'no'},'od',0:3))(3)
nb =
6
% Obtain up-down info
>updown = flgeomud(b)
updown =
1 1 1 0 0 0
0 0 0 1 1 1
% Print normal vector information as obtained by flgeomfd and flgeomfn commands
>for i = 1:nb
[coord,facederiv{i},c] = flgeomfd(b,i,[0.5;0.5]); % Obtain face derivative at mid-surface
normal{i} = flgeomfn(facederiv{i}); % Obtain normal from face-dreivative
end
>normal
normal =
[0 0 1] [1 0 0] [0 1 0] [0 0 1] [1 0 0] [0 1 0]
% The above information lacks direction. Combine with up-down information and
% convert all inward pointing normals into outward pointing normal vectors
>for i = 1:nb
if (updown(1,i)==1) % if inward pointing normal
normal{i}=-normal{i}; % make it outward pointing normal
end
nx(i) = normal{i}(1); % Store in separate vectors (optional)
ny(i) = normal{i}(2);
nz(i) = normal{i}(3);
end
>normal
normal =
[0 0 -1] [-1 0 0] [0 -1 0] [0 0 1] [1 0 0] [0 1 0]
=============
Now that all the outward normal vectors of all the faces have been attained, the total heat energy escaping the box through walls can be easily calculated.
Hope this helps.
Please login with a confirmed email address before reporting spam
Thanks
Jyani
Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.
Suggested Content
- BLOG Understanding Higher-Order Diffraction
- KNOWLEDGE BASE Error: Failed to compute elastoplastic strain variables
- KNOWLEDGE BASE Two-Phase Flow Modeling Guidelines
- BLOG Modeling Electromagnetic Waves and Periodic Structures
- FORUM Accessing Stiffness Matrices and Internal Force Vectors for Incremental Solution in COMSOL