Trouble with finding appropriate conditions for logic

I have a blunt body (Apollo capsule) modeled. I wrote a function to read in the vertices and the normals from the STL of this model. Given a velocity vector, I want to be able to figure out which facets of the capsule are in the shadow region (region where there is relatively no flow) and also any facets not on the heat shield (the facets that are not part of the cone and the radius at the top). I tried writing conditions down using the angle between the velocity vector and the normal vector of each facet. Right now I'm coloring these facets red so I can see if the conditions are detecting the right facets. The problem that I'm having is that some of the facets of interest are getting detected and some areas are not.

Included are some images of what the capsule looks like from the side and in 3d space and also the problem I am running into. The red line indicates the velocity vector pointing up. The whole underside of the capsule should be selected as well as the lower half of the rounded edge (the edge between the large heat shield and the cone sides).

[Fredrick_jong]

Hi,

I was wondering, how did you manage to plot the Apollo capsule on MATLAB? I am currently trying to do that but I am unable to. Your assistance will be greatly appreciated. Many thanks in advance.

Answer 1

[Cedric]

Could you have an issue with the orientation of your facets, e.g. by coding vertice IDs clockwise when they should be coded counterclockwise or vice versa? That would make part of your normals point in the wrong direction which would change the sign of the scalar product with the velocity vector (if it's what you use to detect facets in the shadow region).

Without seeing your code, there isn't much more that comes to my mind.

EDIT: here are a few hints after looking at your code.

• The red line that you display vertically is different from the "free-stream vector" b that you define in facet_sort.m.
• You rotate vertices in Apollo.m but you forget to rotate normals .

These two combined explain the angle. If you set b=[0 0 1] in facet_sort.m for example, you will see that the heat shield will be displayed in red.

Several for loops could be eliminated if you had efficiency requirements. E.g. in facet_sort.m, you could replace the computation of theta (the whole for loop) by:

 flagsShadow = normals * b.' < 0 ;   % Or > 0 depending the orientation
                                     % your facets.

and propagate these flags back into PlotFacets.m to select the color, instead of performing tricky computations with angles.

Hope it helps!

Cedric

[Harold]

Unfortunately I have no control over the IDs of the vertices as these are read from the STL file. Below is the condition that I am using for the test. Here is a link to the code that I'm working with. I've included a read me file for further explanation.

Note: (k+2)/3 is used to calculate the numbers 1,2,3... from k (which goes starts from 1 and steps every 3, i.e. 1, 4, 7,...).

http://www.filehosting.org/file/details/412263/Apollo_Capsule.rar

 if pi/2 < theta((k+2)/3) < 1
   color = 'r';
 else
   color = 'w';
 end

[Cedric]

Please have a look at the edit.

[Harold]

Wow, thank you. Can't believe I made that mistake with not rotating the normals and also the incorrect velocity vector. I made the suggested corrections and everything looks good except for one issue that I can't seem to figure out. Some facets on the underside of the cone are colored red. This is not incorrect in terms of code execution. However for the calculations that are to be performed on this vehicle only facets that are part of the heat shield (big spherical-like portion and the rounded edge) are only considered. I would probably have to add in another condition in facet_sort.m but I am unsure of what else I could use that I already know from the data.

[Cedric]

When you load your data, the z axis is an axis of symmetry and the heat shield is facing downwards, At this moment (before rotation), you can flag facets from the shield or from anywhere else based on the z component; then you can intersect these flags with shadow flags to target more specifically your facets.

You could for example add the following code just before alpha=-60 in Apollo.m:

 hBoundary = 0.6 ;  % Adapt it to your needs.
 flagsShield = any( reshape( vertices(:,3), 3, [] ) < hBoundary ) ;

Here, flagsShield is a vector of logical whose length equals the number of facets and whose value for each facet is true if all vertice z componenent are below 0.6 (just after loading the STL, which seems to target the shield.. but you will want to adjust this value). You then just pass flagsShield to PlotFacets() and use it to target/detect relevant facets.

Cheers,

Cedric

PS: please mark the question as answered if it is. I will go on answering/commenting if there is activity anyway.

[Harold]

Thank you. Clever trick using the z values before rotation. Now I just need to work on the intersection. I'll reply back once I finish.

[Cedric]

if you use vectors of logical for flagging facets, then you just want to apply logical operators

 flagsInRed = ~flagsShadow & flagsShield ;

would be a vector flagging facets that are both in the shield and not in the shadow.

[Harold]

That worked to some degree. The cone on the top is getting selected for some reason even though the z coordinates are greater than 0.65 (the value that I'm testing against). Another problem is that for greater attack angles, the sides are getting picked up. This is odd since the intersection should of taken care of this.

0 attack angle

-30 attack angle

-60 attack angle

Answer 2

[Cedric]

Which method did you chose finally to..

• Flag facets in the shadow.
• Flag facets from the shield.
• Define colors.

Could you paste the code? Here it seems that you have something like a union of the upper part and facets exposed to flow.

[Harold]

Sorry about that. I decided to define red facets as the facets that I want to keep for future calculations such as pressure, temperature, etc. I used both methods flags for the facets which are then passed into PlotFacets.

Apollo.m

 hBoundary = 0.65 ;  % Adapt it to your needs.
 flags_Shield = any( reshape( vertices(:,3), 3, [] ) < hBoundary );

 alpha = 0;        %test value to check and see if facet selection works

 if alpha ~= 0
 [vertices] = rot_trans(vertices(:,1), vertices(:,2), vertices(:,3), alpha);
 [normals] = rot_trans(normals(:,1), normals(:,2), normals(:,3), alpha);
end

 [flags_Shadow] = facet_sort(facet_num, vertices, normals);

 PlotFacets (vertices(:,1), vertices(:,2), vertices(:,3), flags_Shadow,  flags_Shield', handles)

PlotFacets.m

 figure(2)
 axis vis3d
 vcnt = length (vx) - 1;

 %facets that are not going  to be used in future calculations
 flagsInWhite = ~flags_Shadow & flags_Shield; 

 for k = 1: 3: vcnt
    if flagsInWhite((k+2)/3) == 1
        color = 'w';
    else
        color = 'r';
    end
  fill3 ([vx(k), vx(k+1), vx(k+2), vx(k)], [vy(k), vy(k+1), vy(k+2), vy(k)], [vz(k), vz(k+1), vz(k+2), vz(k)], color);
  if k == 1
        hold;
  end
end

 if k > 1
   axis equal;
   hold;
end

 ylabel('y')
 xlabel('x')
 zlabel('z')

[Cedric]

Ok, so if you want to flag in red your region of study which is, up to what I understood, the region of the shield that is not in the "shadow region", then you want to define

 flagsInRed = ~flags_Shadow & flags_Shield ;

and change your code in PlotFacets into

 if flagsInRed((k+2)/3)
   color = 'r';
 else
   color = 'w';
 end

It seems irrelevant just to change a variable name, but actually defining flagsInWhite would mean taking the negation of (~flags_Shadow & flags_Shield) which is (flags_Shadow | ~flags_Shield) ;-)

Now this won't work unless you set the boundary back to 0.6. The reason is that 0.65 captures one or two vertices of these long triangles that constitute the cone, and take the cone as being part of the shield. You could also change the logic behind the detection of polygons from the shield, e.g. by using all() instead of any():

 flags_Shield = all( reshape( vertices(:,3), 3, [] ) < hBoundary );

which would require all vertices to be below the boundary for a facet to be flagged as being part of the shield.

As a side note, I would also change some loops in your code so they use a facetID rather than a vertexID. This would avoid computing the facet ID using (k+2)/3 which is a positive integer because "MATLAB is nice enough to understand that it is likely to be one", but which would be a float in other languages. Starting from a facetID and computing the IDs of vertices by multiplication surely leads to positive integer IDs. But this is more a question of preference I guess..

[Harold]

I changed the variable name back to flagsInRed and changed the code in PlotFacets.m to what you have in your last post. I kept any() for the logic detection as all() did not give desired results. I've included images below to show what the code is giving and what I am looking for. I've should of done this from the get go but didn't think about it. To clear any confusions, white is defined as facets that I want to use for calculations.

[Cedric]

Ah ok, then be sure to have hBoundary=0.6 and not 0.65, and it is just about permuting r/w:

 flagsInWhite = ~flags_Shadow & flags_Shield ;

and

 if flagsInWhite((k+2)/3)
   color = 'w';
 else
   color = 'r';
 end

I made these modifications really fast on the code that you posted yesterday, and it's working.

[Harold]

Works perfect. Thank you very much. I'll be sure to include you in my acknowledgments when I finish this program.

[Cedric]

My pleasure, I'm glad it works!

[Harold]

Cedric, I am working on another bit of code that uses the logic indexes from Flow flags. To help me troubleshoot my code I made some changes. In the facet_sort function I renamed flagsShadow to

Flow_flags = normals * b.' < 0 ;

which returns the logic indexes of all the facets being hit by the oncoming flow.

I also renamed flagsInWhite to

ShieldFlow_flags = Flow_flags & Shield_flags';

This finds the intersection of the Flow_flags and the Shield_flags. This intersection should then return the flags of the facets that are only in the flow that are part of the shield. I then use ShieldFlow_flags in the if/else statement of the plot_facet function to fill the facets with a bluish color. This works perfectly.

While testing code that I developed to calculate the intersection points (points in red; hard to see so download and zoom in) of the facets given a yz plane located at some x value, I discovered something weird. As of right now, my intersection function calculates all the intersection points of every facet. As shown below

Since I am only interested in the intersection points of the facets in the flow and are the shield, I passed vertices(ShieldFlow_flags,:) to my intersection function instead of vertices. On the plot, I noticed that some intersection points were missing.

So I decided to pass vertices(ShieldFlow_flags,:) into plot_facets. I found that some facets were not getting filled in.

So the question is this, why is it if I use the entire vertices matrix and fill in facets with ShieldFlow_flags everything works but if I use vertices(ShieldFlow_flags,:) there are some facets not getting filled? I have been trying to figure this out for the past few days but can not understand what the problem is. I understand if you're busy. This project is something for myself, so there is no time constraint.

[Cedric]

Well, without running your code, my first guess is that you are indexing your array of vertices whose size is 3*nFacets x 3 (= nVertices x 3 ) with a logical index for facets whose size is nFacets.

[Harold]

I checked the size of the logic indexes and the vertices and you were right. The size of the logic indexes is exactly 3 times smaller than that of the vertices. This is because Flow_flags is calculated using normals. The size of normals will always be 3 times smaller than the vertices since each row or normals describes the location of the normal. So what I did was output Flow_flags from face_sort and used the kron function to replicate every row of Flow_flags three times so that I get the same size as the vertices matrix. This worked to some extent. I have a few intersection points that appear on the small rounded cone part at the top of the capsule. This doesn't do this when alpha is set to 0 though. I'm not sure why this happens but I know that I am on the right track. By the way, the variable xsym is the x value of where the plane of symmetry is for the model. You can change this to change the x location of the yz cutting plane. Just run Apollo from the editor window.

http://www.filehosting.org/file/details/415545/Apollo.rar

I just realized that if the rounded cone comes past the heat shield, there are flags that get picked up. If you set view([-90 0]) and change the alpha angle you can see this happen.

Answer 3

[Harold]

I believe I have just figured out what I was doing wrong. I should of been using the replicated ShieldFlow_flags for the row indexes of vertices. This works now. I now only get red dots on the heat shield facets that are exposed to the flow.

[Cedric]

The key is: if it is not obvious whether a variable/function applies to vertices or to facets, then you should name this variable/function so it explicitly gives this information. Here are a few hints..

• For loops: name a loop index fId if it applies to facets and vId if it applies to vertices, instead of using i, j, k, etc.
• Build functions for converting vectors of indices from vertices to facets and vice versa. You could even use short names, e.g. f2v() and v2f():

 function vId = f2v(fId)
    % Return a nx3 array of vertices IDs, where n is the number of elements
    % of fId (vector of facet IDs).
    vId = [3*fId(:),3*fId(:)+1, 3*fId(:)+2] ;
 end
 function fId = v2f(vId)
    % ...
    ...
 end

• Think about "vector" ways to process your data. For example, If you had to build a flag enabling/disabling all facets that contain vertices from a list of vertices IDs vId_list, it would be much more efficient/clear to perform the following (using functions defined above):

 fId_list = v2f(vId_list) ;
 fEnabled = true(nFacets, 1) ;
 fEnabled(fId_list) = false ;

than to build a complicated for loop using i as index, rounding i/3, etc.

[Harold]

Thanks for the tips. I will definitely have to incorporate these.

[Cedric]

In fact, you could have the transpose of what I built above for vId, so you can exploit linear indexing:

 function vId = f2v(fId)
    % Return a 3xn array of vertices IDs, where n is the number of elements
    % of fId (vector of facet IDs).
    vId = [3*fId(:),3*fId(:)+1, 3*fId(:)+2].' ;
 end

I don't know if you are familiar with linear indexing.. if not, look at the following:

 >> A = [1 2; 3 4] ;
 >> A(:)
 ans = 
        1
        3
        2
        4

Arrays are saved in memory in "column". A(:) means get all elements of A following the linear order in memory, which means that you get the entire column 1 first and the the entire column 2. This means for example that is v is a vector, v(:) is always the column version of this vector, even if v is a row vector.

So having vId as a 3xn matrix means that if you access it linearly with vId(:) you get a column vector with the 3 vertices of facet 1 of fId first, then the 3 vertices of facet 2 of fId, etc .. which can be useful.