Recently, I created a project called ForCAD, which generates Bezier, Rational Bezier, B-spline, and NURBS curves, surfaces, and volumes. ForCAD functions as a mesh generation tool, applicable for finite element method (FEM) and isogeometric analysis (IGA).

This looks promising. I would suggest you add the following in future releases (particularly if you want to support IGA).

- Derivatives of objects (curves, surfaces, volumes) along with derivatives of basis functions
- Parametrization of objects (knot vectors)
- Merging of curves and surfaces with at least C1 continuity.
- Extraction of piecewise Bezier objects from NURBS
- degree elevation and deflation (knot insertion).

There are a host of other things that will be required (and probably several years work) if you want to match openCascade or the SINTEF GO tools but this is a good start.

Thanks @rwmsu for your feedback.

I have implementations for derivatives; this will be added soon.

I tried achieving C^0 continuity for multiple patches, but achieving higher continuity for multiple patches seems not easy. Actually, I’ve been looking for such an algorithm for a while. If you know of any algorithm or references for that, I would be happy to learn about it.

Yes, this is very much needed. Degree elevation for Bezier and rational Bezier curves is supported. It must be extended for surfaces and volumes. Degree elevation and reduction, as well as knot insertion and removal, are on the to-do list.

Do you mean a derived type for knot vectors? Actually, I have a more complicated implementation for NURBS, which sometimes makes development very difficult. I thought of simplifying it for ForCAD. But a derived type for knot vectors is a good idea! I will add it to the to-do list.

This is also a good idea. Since I haven’t used it within IGA, I never thought to implement it. I will add it to the to-do list.

Yes, there are many features that can be added to ForCAD. If there is any interest in contribution, I would be happy.

By parameterization, I mean generating the normalized (or parameterized) knot vector based on some distribution in physical space, normally some measure of arc length for curves. The idea is to find something approaching an optimum knot vector spacing that minimizes errors in the resulting object. For surfaces and volumes it gets more complicated and will involve techniques like solving a harmonic (elliptic) system of equations to get something near to an optimum distribution. See the parameterization code based on the work of Michael Floater in the SINTEF GOtools software. I do recommend a separate knotVector class. About 15 years ago I started writing my own Fortran NURBS package to support IGA but gave it up when I couldn’t convince anybody to fund me to complete it. I’m retired now so maybe I’ll give it another look

This is great, the Fortran open source community keeps growing! I’ve been looking on and off for a package to calculate FEA element volumes for a while now, seem like this would be able to do that?

Thank you, @rwmsu. It’s great to hear that you’ve worked on NURBS and IGA. I also have an implementation for IGA, but at the moment, I don’t have time to create a repository for it. I will check out SINTEF GOtools.

@Chuckyvt , ForCAD is currently focused on describing geometries using NURBS. However, if you have FEM or IGA software, it can be used as a preprocessor to generate meshes. ForCAD provides you with element connectivity.

ForCAD v0.2.0 now supports `knot insertion`

, `degree elevation`

, and `calculation of derivatives`

for NURBS curves, surfaces and volumes. To help test the functionalities and improve the library, please submit any NURBS examples if you already have as pull requests (PRs). Your contributions will greatly enhance ForCAD. Thank you!

ForCAD now supports knot removal, embedding meshes into NURBS, retrieving IGA element connectivity, and includes some predefined NURBS shapes such as circle, half circle, tetragon, hexahedron, 2D ring, half 2D ring, 3D ring, half 3D ring. It also includes rotatation and translation of geometry and control geometry, along with visualization using PyVista. Below is an example of embedding meshes into a NURBS volume and visualizing it using PyVista.

A logo has also been designed using ForCAD and visualized using ParaView.

Additionally, I have added new examples that convert NURBS surfaces to PPM images using ForImage and colorize them using ForColormap.

@Ali

That’s great news!

I am the developer of some Fortran codes that make use of the SISL and forSISL (thank you @rwmsu ) libraries to model and mesh Turbomachinery and Aerospace components for CFD applications…

Given that SISL is no longer maintained ForCAD seems perfect for a future replacement!

A capability to project x,y,z points on a surface (i.e. getting u and v parameters that are close enough to the initial point ) would be particularly useful!

This is an operation that is carried out with high frequency in CFD meshing, e.g. connecting a wing to the aircraft body or a compressor blade to the hub/root platform…

SISL has got that in two flavours: a high accuracy routine and a very fast one… the latter being particularly useful and time-saving if you have got thousands of points to project!

Thank you for sharing your project code!

Agree totally. It would be very nice if the (very) long range goals of forCAD would be to match the functionality of openCascade or the SINTEF GoTools. I say very long range because to do the job right will take a lot of time and probably more than one person envolved in the development.

Edit. I think a very real feature that needs to be addressed from the start are IGES and STEP readers and writers. These do not have to support the full IGES or STEP standard but should include NURBS/Bezier support or a subset such as the NASA IGES subset

Also a version of the full IGES 5.3 standard from 1996 is available (I presume legally) from the following Lawrence Berkeley labs web site.

Yeah, I would say that this would be an amazing library to have in Fortran although I suspect getting feature-parity with something like OpenCASCADE is very hard to achieve.

A subset of STEP and IGES is probably easier to achieve, OpenCASCADE has a relatively complete implementation which could be used as inspiration. Reading the STEP ISO Standards is somewhat painful not to mention pricey!

@Rob777, I’m glad you find ForCAD interesting, and I’ll be even happier if you use ForCAD for real simulations.

I’ve added a method `nearest_point()`

for `nurbs_surface`

. It’s now available in the dev branch. Additionally, I’ve created an example: nearest_point_2d.f90. I’m looking forward to your feedback to ensure it meets your requirements!

@rwmsu Thanks for sharing the information! I have a to-do list in ROADMAP.md where I’m gathering ideas. IGES and STEP files are indeed on the list (I think), and I believe this addition is necessary. Currently, I only use ASCII vtk files to export generated objects.

Initially, I created ForCAD for use within FEM and IGA simulations. I believe ForCAD is now capable of being used within FEM and IGA codes. If I find the time, I will add an example for FEM and IGA.

Since I develop ForCAD and other Fortran codes in my free time, adding many features in a short time isn’t easy, as you mentioned. It would require many copies of me;)!

I’ll be happy for any contributions and your feedback!

@Ali thanks for the example code! It looks nice just by reading it even if I still have to find the time to download it and actually try it…

What is exactly the id variable for?

I would expect nearest to contain the u,v coordinate one is looking for …

To clarify, I use the following notation:

- \mathbf{X}_t \in \mathbb{R}^{n_t \times 2} represents the coordinates in parametric space, which typically range from 0 to 1.
- These points are mapped from parametric space to physical space, denoted as \mathbf{X}_g \in \mathbb{R}^{n_g=n_t \times 3}.

The `nearest_point(given_point, nearest, id)`

method identifies the nearest point to `given_point`

in physical space, which refers to the nearest point on the surface. `id`

is the index of \mathbf{X}_g, indicating which point on the surface is nearest to the `given_point`

. If u and v are the coordinates of the parametric space \mathbf{X}_t, I can easily extend the method to provide these coordinates.

Thank you for your kind explanation but still not that clear to me …

You haven’t specified if `given_point`

is Xt or Xg… the same applies to `nearest`

Let’s keep it simple …

I have got a surface S (it could be the upper side of a wing) embedded in the 3-dimensional space x,y,z…

I also have a point P belonging to x,y,z located close enough to this surface…

What I’d like to know are the parametric coordinates u,v of the point p0 that lies on S and is the closest possible to P…

How should I use the routine nearest point?

Something like this

```
real(wp), dimension(3) :: P
real(wp), dimension(2) :: p0
integer :: id
P=[whatever_x, whatever_y,whatever_z]
call S%nearest_point(P,p0,id) ! what is id?
```

Thank you so much

I have updated the `nearest_point()`

method and the corresponding example:

```
!> Find the nearest point on the surface to a given point
! nearest_Xg: Coordinates of the nearest point on the surface (optional)
! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional)
! id: id of the nearest point (optional)
call shape%nearest_point([2.0_rk, 3.0_rk, 5.0_rk], nearest_Xg, nearest_Xt, id)
print *, 'Nearest point on the surface:', nearest_Xg, 'with parametric coordinates:', nearest_Xt, 'and id:', id
```

Output of the example:

`Nearest point on the surface: 2.0 3.0 0.0 with parametric coordinates: 1.0 1.0 and index: 900`

This indicates that the 900th point on surface S has the coordinates (x=2, y=3, z=0), which are the nearest to the point P=(x=2, y=3, z=5). The corresponding parametric coordinates of (x=2, y=3, z=0) are (u=1, v=1).

(All points describing the surface in physical space are stored in a matrix X_g, and the corresponding coordinates of these points in parametric space are stored in a matrix X_t. The id indicates which row of these matrices.)

I hope I’ve clarified it, but if it’s still unclear, I can provide a visual example with more details.

Edited:

Regarding performance, I’ve improved the `nearest_point()`

subroutine with OpenMP and performed a simple test (nearest_point_2d_bench.f90). Here are the results on my system with `gfortran 14.1.0`

:

```
fpm run --example nearest_point_2d_bench --profile release --flag "-ftree-parallelize-loops=16"
```

Dear @Ali,

I fear I have not fully understood how the nearest_point algorithm operates.

It seems that, given a surface S and a 3D point P (possibly located very close to that surface but not lying on it), the routine will find the closest point to P choosing from a set of point…

What this set of point should be is not clear to me…

The way I’d expect this routine to perform is to find the closest point to P with the only constraints that the solution point lies on the surface and possibly with an arbitrary precision…

Why are you limiting the search to a set of points to choose from?

Any point on the surface could be a valid solution…

Please for a reference have a look at functions **s1958** and **s1954** in the SISL library…

Thank you for your time and work