[Mars409]

Thanks.

Still in iteration of inspecting codes and correcting inconsistencies -- they crop up naturally as concepts evolve. One question: which class.method calculates the point coordinates to be printed out for spline edge defined by OnCurve class? OnCurve.representation() provides the keyword "spline" by default but I cannot find the codes to output point coordinates. OnCurve.discretize(self, param_from, param_to) seems the ideal candidate, but where is it called when writing the blockMeshDict?

Oh wait, is it CurveEdgeBase.description() https://github.com/damogranlabs/clas...s/curve.py#L33

mesh.py:            output.write(self.edge_list.description)

class EdgeList:
    """Handling of the 'edges' part of blockMeshDict"""
    ....
    @property
    def description(self) -> str:
        """Outputs a list of edges to be inserted into blockMeshDict"""
        out = "edges\n(\n"

        for edge in self.edges:
            out += edge.description + "\n"

        out += ");\n\n"

        return out

@dataclasses.dataclass
class Edge(ElementBase):
    """Common stuff for all edge objects"""

    vertex_1: Vertex
    vertex_2: Vertex
    data: EdgeData
    @property
    @abc.abstractmethod
    def description(self) -> str:
        """string description of the edge to be put in blockMeshDict"""
        # subclasses continue from here
        return f"\t{self.representation} {self.vertex_1.index} {self.vertex_2.index} "

    @property
    def representation(self) -> EdgeKindType:
        """A string that goes into blockMesh"""
        return self.data.kind

class EdgeData(ElementBase):
    """Common operations on classes for edge creation"""

    kind: EdgeKindType

class OnCurve(EdgeData):
    """An edge, snapped to a parametric curve"""

    kind: EdgeKindType = "curve"

    def __init__(self, curve: CurveBase, n_points: int = 10, 
                         representation: EdgeKindType = "spline"):
    ....
        self._repr: EdgeKindType = representation

    @property
    def representation(self) -> EdgeKindType:
        return self._repr

 class CurveEdgeBase(Edge, abc.ABC):
    """Base class for edges of any curved shape,
    defined as a list of points"""

   @property
    def description(self):
        point_list = " ".join([vector_format(p) for p in self.point_array])
        return super().description + "(" + point_list + ")"

along with OnCurveEdge.point_array()? https://github.com/damogranlabs/clas.../curve.py#L83?

@dataclasses.dataclass
class OnCurveEdge(CurveEdgeBase):
    """Spline edge, defined by a parametric curve"""

    data: edges.OnCurve

...
    @property
    def param_start(self) -> float:
        """Parameter of given curve at vertex 1"""
        return self.data.curve.get_closest_param(self.vertex_1.position)

    @property
    def param_end(self) -> float:
        """Parameter of given curve at vertex 2"""
        return self.data.curve.get_closest_param(self.vertex_2.position)

    @property
    def point_array(self) -> NPPointListType:
        return self.data.discretize(self.param_start, self.param_end)[1:-1]

class OnCurve(EdgeData):
    """An edge, snapped to a parametric curve"""

    kind: EdgeKindType = "curve"
   ...
    def discretize(self, param_from: float, param_to: float) -> NPPointListType:
        return self.curve.discretize(param_from, param_to, self.n_points + 2)

https://github.com/damogranlabs/clas.../edges.py#L163

[kandelabr]

Thou Art Not Mistaken! It's in the CurveEdgeBase.


I am curious though (and a bit worried), how come you're inspecting that code so meticulously and into such gory details?

[Mars409]

Because I am not sure how exactly to define the edges for the sketches after wrapping onto the cylinders. Until two days ago, I was not paying attention to this aspect, and only focused on wrapping the points onto cylinder and projecting clamp lines onto the cylinders. These I got them worked out completely, finally, and self-consistently. But then I realize the point-to-point lines in the sketches cannot be left as default straight edges but instead every one of them must become edges either as arc or spline.

Now, I finally understand how the spline edge's list of points will be generated and n_points is the number of points inserted into the curve between the pair of ADJACENT vertices, and NOT between the two FAR ENDS of the entire curve -- this was THE KEY QUESTION I need answer for, and now I have found it after going thru the codes. So, in a nutshell, when an edge is defined using an On_Curve which has been initialized with an n_point value, there will be exactly n_point+2 points in the spline list of points in BlockMeshDict. The n_point+2 points are NOT spread all over the enter curve used to initialize the On_Curve object.

Finally, I suggest a short note goes into the airfoil tutorial to clarify what n_points is. It can save everyone else that comes after this a lot of time.

Even now, I am not 100% comfortable with defining edges within the sketches by splines b/c I don't have very good control over the error of departure from the cylinders. I am next going to look into possibly using arc edges instead.

[kandelabr]

Originally Posted by Mars409

Finally, I suggest a short note goes into the airfoil tutorial to clarify what n_points is. It can save everyone else that comes after this a lot of time.

Noted! Will do.

Originally Posted by Mars409

Even now, I am not 100% comfortable with defining edges within the sketches by splines b/c I don't have very good control over the error of departure from the cylinders. I am next going to look into possibly using arc edges instead.

Try the Origin edge - it will position the arc to pass through both points and have a center on the origin. For small angles, the error shouldn't be too big even if it should actually be an ellipse/spline/whatever...

[Mars409]

It turns out the OF Foundation's BlockMesh angle-and-axis specification for edge fully supports specifying a straight line wrapped onto a cylinder, aka helix. Their user guide is silent on this, but their C++ codes shows:

https://github.com/OpenFOAM/OpenFOAM.../arcEdge.C#L75

void Foam::blockEdges::arcEdge::calc(const scalar theta, const vector& axis)
{
    ...
    const vector dp = p1_ - p0_;

    const vector pM = (p0_ + p1_)/2;
    const vector rM = normalised(dp ^ axis);

    const scalar l = dp & axis;

    const vector chord = dp - l*axis;
    const scalar magChord = mag(chord);

    centre_ = pM - l*axis/2 - rM*magChord/2/tan(theta/2);
    axis_ = axis;
    theta_ = theta;
    length_ = l;
}

and

https://github.com/OpenFOAM/OpenFOAM...arcEdge.C#L133

Foam::blockEdges::arcEdge::arcEdge
(
    const dictionary& dict,
    const label index,
    const searchableSurfaces& geometry,
    const pointField& points,
    Istream& is
)
:
    blockEdge(dict, index, points, is),
    p0_(points_[start_]),
    p1_(points_[end_]),
    centre_(),
    axis_(),
    theta_(),
    length_()
{
    token firstToken(is);
    is.putBack(firstToken);

    if (firstToken == token::BEGIN_LIST)
    {
        ...
    }
    else
    {
        const scalar theta = degToRad(readScalar(is));
        const vector axis(is);
        calc(theta, normalised(axis));
    }
}

and

https://github.com/OpenFOAM/OpenFOAM...arcEdge.C#L152

Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
{
    ...

    const scalar t = theta_*lambda;
    const vector r1 = p0_ - centre_, r2 = axis_ ^ r1;

    return centre_ + r1*cos(t) + r2*sin(t) + axis_*length_*lambda;
}

There's only one problem: the length function has a bug:

Foam::scalar Foam::blockEdges::arcEdge::length() const
{
    const vector r1 = p0_ - centre_;

    // Length of a helical segment
    return theta_*sqrt(magSqr(r1) + sqr(length_));   
    #             ^^^^^^^^^^^^^^^ BUG ABOVE. 
    # Correction: sqrt(sqr(theta_)*magSqr(r1) + sqr(length_))
}

The previous version, v11, has even worse bug:

Foam::scalar Foam::blockEdges::arcEdge::length() const
{
    const vector r1 = p0_ - centre_;

    // Length of a helical segment
    return degToRad(theta_*sqrt(magSqr(r1) + sqr(length_)));  
    #      ^^^^^^^^^^ WORSE BUG
}

[Mars409]

My next question is: Does the edge need to be updated in each iteration of optimization? My guess is Yes, b/c the internal angles of the face and the normal to the face do depend on the edge.

The SPLINE edge defined thru OnCurve has no problem with updating without needing user code, b/c the angles at any junction along a curve that constructs the OnCurve can be computed within the cb codes. (I don't know yet whether this is what actually happens ... )

The ARC edge, otoh, seems to need user code to adjust the angle and axis for every iteration of optimization if it requires information as to the edge's geometry.

So I have to investigate by TRACING thru the cb codes AGAIN ... But maybe you can tell the answers me straight away.

=====
ANSWER FOUND : No, b/c the internal angles of the face and the normals to the face triangles depend only on the points that form corners of the sides of the hex, as if the edges were perfectly straight lines, and NOT on the actual edges, however bent they may actually be. So, the edge definitions are completely irrelevant to the optimization, and can be specified just once for all edges.
https://github.com/damogranlabs/clas...ze/cell.py#153

    def get_edge_lengths(self) -> FloatListType:
        points = self.points

        return np.array([f.norm(points[edge[1]] - points[edge[0]]) for edge in self.side_indexes])

    @property
    def quality(self):
    ...
                    c2c = center - neighbour.center

                c2cn = c2c / np.linalg.norm(c2c)

                angles = 180 * np.arccos(np.dot(self.get_side_normals(i), c2cn)) / np.pi
                quality += np.sum(q_scale(1.25, 0.35, 0.8, angles))
                ### cell inner angles
                quality += np.sum(q_scale(1.5, 0.25, 0.15, abs(self.get_inner_angles(i))))

            ### aspect ratio: one number for the whole cell (not per side)
            edge_lengths = self.get_edge_lengths()
            side_max = max(edge_lengths)
            side_min = min(edge_lengths) + VSMALL
            aspect_factor = np.log10(side_max / side_min)

            quality += np.sum(q_scale(3, 2.5, 3, aspect_factor))

Even the edge_lengths actually are not edge lengths but straight line distances between connected points.

Below codes from HexCell shows methods get_side_normals() and get_inner_angles() uses point coordinates only.
https://github.com/damogranlabs/clas...e/cell.py#L210

class HexCell(CellBase):
    """A block, treated as a single cell;
    its quality metrics can then be transcribed directly
    from checkMesh."""
...
    def get_side_normals(self, i: int):
        side_center = self.get_side_center(i)
        side_points = self.get_side_points(i)

        side_1 = side_points - side_center
        side_2 = np.roll(side_points, -1, axis=0) - side_center

        side_normals = np.cross(side_1, side_2)

        nnorms = np.linalg.norm(side_normals, axis=1) + VSMALL
        return side_normals / nnorms[:, np.newaxis]
...
    def get_inner_angles(self, i: int):
        side_points = self.get_side_points(i)

        sides_1 = np.roll(side_points, -1, axis=0) - side_points
        side_1_norms = np.linalg.norm(sides_1, axis=1) + VSMALL
        sides_1 = sides_1 / side_1_norms[:, np.newaxis]

        sides_2 = np.roll(side_points, 1, axis=0) - side_points
        side_2_norms = np.linalg.norm(sides_2, axis=1) + VSMALL
        sides_2 = sides_2 / side_2_norms[:, np.newaxis]

        angles = np.sum(sides_1 * sides_2, axis=1)
        return 180 * np.arccos(angles) / np.pi - 90

======

NEVERTHELESS, since the ARC's angle specification depends on the two endpoints' coordinates, which are changed by optimization, their face.add_edge(corner, edge_data) calls DO need to be performed AFTER optimization.

[kandelabr]

When I was 'transcribing' foundation's angle-and-axis code there was no mention of cylinder geometry or whatever. It is just a different representation of an arc. When origin and the two points don't form a circle (two different radii, for instance), origin is moved so that a proper arc can be formed. This is what classy_blocks does at the moment - wrapping an arc onto a cylinder would need an ellipse, though - not sure if OF has that at the moment.


You're right, edges have to be redefined after optimization. The latter doesn't take into account edges anyway (just blocks with straight lines). When arc edge's vertices are moved, edge definition stays the same - move a lot and it will break your mesh.


If you move a spline edge's vertex along that spline, that will be fine. When you move away from it, it probably won't work well.

[Mars409]

Originally Posted by kandelabr

When I was 'transcribing' foundation's angle-and-axis code there was no mention of cylinder geometry or whatever. It is just a different representation of an arc. When origin and the two points don't form a circle (two different radii, for instance), origin is moved so that a proper arc can be formed. This is what classy_blocks does at the moment - wrapping an arc onto a cylinder would need an ellipse, though - not sure if OF has that at the moment..

The foundation's arcEdge.C code snippets I showed above show that the foundation's angle-and-axis ARC edge takes care of wrapping the line onto a cylinder. The l (i.e. the letter after k, not numeral one, though the font may be indistinguishable from the numeral 1) is the component of displacement along the cylinder axis. Where l=0, it's an arc along a circle. Where l > 0, it's a helix. The header file arcEdge.H explicitly says this. https://github.com/OpenFOAM/OpenFOAM.../arcEdge.H#L32

    If the angle-axis specification is used, the axis may not be not
    perpendicular to the line connecting the end points. In that case, the axis
    is considered to be that of a cylinder, and the arc represents a portion of
    a helix on the surface of that cylinder.

and https://github.com/OpenFOAM/OpenFOAM...e/arcEdge.H#79 :

        //- The axial length of the arc/helix
        scalar length_;

[kandelabr]

Yes but it is still just a circular arc. The 'geometry' argument is not used in the arc code so it can't do anything but an arc anyway.


For a narrow segment an arc will be a good approximation but wider and more 'inclined arcs would become ellipses.


Classy always uses arc definition but adds the original blockMesh definition as a comment above it. You can uncomment that and see if there's a difference.

[Mars409]

Originally Posted by kandelabr

Yes but it is still just a circular arc. The 'geometry' argument is not used in the arc code so it can't do anything but an arc anyway.

For a narrow segment an arc will be a good approximation but wider and more 'inclined arcs would become ellipses.

See position(lambda) function again reproduced, below. Here, the returned position relative to the centre consists of two parts, (1) r1*cos(t) + r2*sin(t), which describes the trajectory of a circle in a plane perpendicular to the axis and passing thru the centre, and (2) axis_*length_*lambda, which describes the simultaneous drift parallel to the cylindrical axis.

This is precisely a helix if length_ > 0, not an arc. It's a arc only if length_ == 0.

https://github.com/OpenFOAM/OpenFOAM...arcEdge.C#L152

Foam::point Foam::blockEdges::arcEdge::position(const scalar lambda) const
{
    ...

    const scalar t = theta_*lambda;
    const vector r1 = p0_ - centre_, r2 = axis_ ^ r1;

    return centre_ + r1*cos(t) + r2*sin(t) + axis_*length_*lambda;
}