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

Please login with a confirmed email address before reporting spam

Hi,

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

7 Replies Last Post Dec 4, 2014, 10:32 a.m. EST
Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Dec 22, 2010, 1:42 p.m. EST
Hi

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
Hi 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

Posted: 1 decade ago Dec 22, 2010, 2:49 p.m. EST
Well, it would be ideally if I could calculate a normal vector in the center of each element pointing outwards. That would be a on the boundary surface. So: All I want to do is calculate a surface vector somewhere on the boundary surface of a 3D geometry that is not given by a parameter but directly by its coordinates.

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.
Well, it would be ideally if I could calculate a normal vector in the center of each element pointing outwards. That would be a on the boundary surface. So: All I want to do is calculate a surface vector somewhere on the boundary surface of a 3D geometry that is not given by a parameter but directly by its coordinates. 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.

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Dec 22, 2010, 3:00 p.m. EST
Hi

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
Hi 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

Posted: 1 decade ago Dec 23, 2010, 3:41 a.m. EST
My geometry can be pretty curvy and pretty long :-)

As you said - (nx,ny,nz) are defined per element. But how can I get this information per element?
My geometry can be pretty curvy and pretty long :-) As you said - (nx,ny,nz) are defined per element. But how can I get this information per element?

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Dec 23, 2010, 4:46 a.m. EST
Hi
try an arrow plot on boundaries with (nx,ny,nz) then in matlab you can retrieve these variables as any others

--
Good luck
Ivar
Hi 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

Posted: 1 decade ago Feb 9, 2012, 1:02 a.m. EST
Calculation of normal vectors on a surface can be tricky. This is how I understand it. Please feel free to correct me.

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.
Calculation of normal vectors on a surface can be tricky. This is how I understand it. Please feel free to correct me. 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

Posted: 9 years ago Dec 4, 2014, 10:32 a.m. EST
These functions seem to be obsolete in the latter COMSOL versions. Any thoughts on how to do this with Matlab Livelink in version 4.3a ?

Thanks
Jyani
These functions seem to be obsolete in the latter COMSOL versions. Any thoughts on how to do this with Matlab Livelink in version 4.3a ? 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.