
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">

<html>

<head>

  <meta http-equiv="Content-Type" content="text/html;charset=ISO-8859-1">

  <title></title>

</head>

<body>

Jeff,<br>

<blockquote type="cite" cite="mid410E2CF6.6060002@cdnorthamerica.com">&nbsp;you

need to intersect the  streamline with the terminating boundary.&nbsp; [...]&nbsp;

Are you sure that your intersection was successful? <br>

</blockquote>

The problem is that for some perfectly legal input the algorithm does not

get there.<br>

<br>

There are two important points in the "while" loop that generates the<br>

points on the stream line:<br>

- at the computation of the next point:<br>

<i>if ((tmp= integrator-&gt;ComputeNextStep(...)) != 0)<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; {<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; /********************************** my code *************************************/<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; if (tmp!=OUT_OF_DOMAIN)<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; break;<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; // compute intersection with edge and try to switch to next

cell<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp; &nbsp; // ....<br>

}</i><br>

<br>

- at the computation of the velocity vector for the newly generated point:<br>

<i>if ( !func-&gt;FunctionValues(point2, velocity) )<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; {<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; retVal = OUT_OF_DOMAIN;<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; memcpy(lastPoint, point2, 3*sizeof(double));<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; break;<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; }</i><br>

<br>

Now for a very simple example the loop exits from the second condition<br>

(I enclose at the end of this message the .vtk file). This means that the

<br>

integration step has taken place succesfully, but the FunctionValues<br>

fails; more exactly, it calls EvaluatePosition, and this fails. Note that<br>

the point is in the current cell - see the attached image (the seed point

is<br>

(0.05, 0.02, 0).<br>

<br>

Each computed point is projected on the current cell. I inserted the following<br>

code right before inserting the point. I wrote a function "Project" because

I <br>

do not really know the meaning of the 0 or -1 returned by EvaluatePosition,<br>

if the point under consideration is outside the cell and its projection onto

the <br>

cell's plane is internal to the cell (or on its boundary).<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp; <br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;<i> if (cell-&gt;GetCellDimension()==2)<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; {<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; double projPoint[3];<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; if (Project(cell, point1, projPoint))<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; {<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; memcpy(point1, projPoint, 3*sizeof(double));<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; }<br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; }</i><br>

&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;<br>

I think that for this case integrator-&gt;ComputeNextStep is successful because<br>

both the initial point and the point+delta/2 fall within the cell; the next

point, <br>

however, gets out of the cell, and FunctionValues somehow detects this.<br>

Which means that I have to hack also around the call of FunctionValues.<br>

<br>

Petru<br>

<br>

-------------------<br>

# vtk DataFile Version 2.0<br>

Sample data for polydata<br>

ASCII<br>

<br>

DATASET POLYDATA<br>

POINTS 5 float<br>

0.0 0.0 0.0<br>

2.0 0.0 0.0<br>

0.0 2.0 0.0<br>

2.0 2.0 1.0<br>

1.5 3.5 1.5<br>

POLYGONS 3 12<br>

3 0 1 2<br>

3 1 2 3<br>

3 2 3 4<br>

<br>

POINT_DATA 5<br>

SCALARS scalars float 1<br>

LOOKUP_TABLE default<br>

0.0<br>

1.0<br>

2.0<br>

3.0<br>

4.0<br>

VECTORS vectors float<br>

1 &nbsp;&nbsp;&nbsp; 1&nbsp; &nbsp;&nbsp;&nbsp; 21<br>

4&nbsp;&nbsp;&nbsp; 2&nbsp; &nbsp;&nbsp;&nbsp; 21<br>

6&nbsp; &nbsp;&nbsp;&nbsp; 3&nbsp; &nbsp;&nbsp;&nbsp; 21<br>

8&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; 4&nbsp; &nbsp;&nbsp;&nbsp; 21<br>

10&nbsp;&nbsp;&nbsp; 5 &nbsp;&nbsp;&nbsp; 21<br>

<br>

</body>

</html>