CAID Python Package

Tutorial

Predefined geometries

square

You can create a square by importing the square function

# ... creates a unit square of degree (2,2) and 4 elements in each direction
from caid.cad_geometry import square
geo = square(n=[3,3], p=[2,2])
geo.plotMesh(MeshResolution=10)
plt.savefig("predefined_geometries_2d_square.png")
# ...

The resulting plot is

_images/predefined_geometries_2d_square.png

circle

You can create a circle by importing the circle function

# ... creates a unit circle of degree (2,2) and (8,4) elements in each direction
from caid.cad_geometry import circle
geo = circle(n=[7,3], p=[2,2])
geo.plotMesh(MeshResolution=10)
plt.savefig("predefined_geometries_2d_circle.png")
# ...

The resulting plot is

_images/predefined_geometries_2d_circle.png

quart_circle

You can create a quarter of circle by importing the quart_circle function

# ... creates a unit quart_circle of degree (2,2) and 4 elements in each direction
from caid.cad_geometry import quart_circle
geo = quart_circle(n=[3,3], p=[2,2])
geo.plotMesh(MeshResolution=10)
plt.savefig("predefined_geometries_2d_quart_circle.png")
# ...

The resulting plot is

_images/predefined_geometries_2d_quart_circle.png

annulus

You can create a annulus by importing the annulus function

# ... creates a unit annulus of degree (2,2) and 4 elements in each direction
from caid.cad_geometry import annulus
geo = annulus(n=[7,7], p=[2,2])
geo.plotMesh(MeshResolution=10)
plt.savefig("predefined_geometries_2d_annulus.png")
# ...

The resulting plot is

_images/predefined_geometries_2d_annulus.png

circle 5 patchs

You can create a 5 patchs circle by importing the circle_5mp function

# ... creates a unit circle with 5 patchs, of degree (2,2) and (4,4) elements in each direction
from caid.cad_geometry import circle_5mp as circle
geo = circle(n=[7,7], p=[2,2])
geo.plotMesh(MeshResolution=10)
plt.savefig("predefined_geometries_2d_circle_5mp.png")
# ...

The resulting plot is

_images/predefined_geometries_2d_circle_5mp.png

Geometry Transformations

The nice thing about Bézier and B-Spline surfaces, is that any affine transformation on the geometry can be done on its control points.

Translation

Let’s start from a unit square and translate it by a displacement vector (1,2). This can be done easily by

# ... creates a unit square of degree (2,2) and 4 elements in each direction
geo = square(n=[3,3], p=[2,2])
# ...

# ... translate the geometry
displ = np.array([1.0, 2.0])
geo.translate(displ)

The result is the following plot

_images/transformations_2d_translate.png

Rotation

Let’s take a unit square and rotate it with and angle of \frac{\pi}{4}

# ... creates a unit square of degree (2,2) and 4 elements in each direction
geo = square(n=[3,3], p=[2,2])
# ...

# ... rotate the geometry
geo.rotate(np.pi/4, axis=2)
# ...

The result is the following plot

_images/transformations_2d_rotate.png

Scaling

In this example, we show how to construct en ellipse starting from a unit circle. We take a unit circle and scale it with (2,4) in each direction

# ... creates a unit circle of degree (2,2) and 4 elements in each direction
geo = circle(n=[3,3], p=[2,2])
# ...

# ... scaling with (2,4) in each direction
geo.scale(2., axis=0)
geo.scale(4., axis=1)
# ...

The result is the following plot

_images/transformations_2d_scale.png

Dive into CAID

CAID is based on a two level objects. The first one is describing a single patch (subdomain) which is the cad_nurbs object for the quadrangular case. The second one (level 0) describes the global geometry as a collection of level one objects: cad_geometry

cad_nurbs

The cad_nurbs object is an extension of the NURBS of IGAKIT but implements some additional informations that are needed for the multi-patchs case, like:

  • orientation: needed for Neumann boundary conditions
  • rational: True if we use the weights. Default value : False

The following script creates a curve, inserts one knot and raises the spline degree by one:

import numpy as np
from caid.cad_geometry import cad_nurbs
from matplotlib import pyplot as plt
from numpy import pi, array

# ... create a 3D curve
def make_curve():
    T = np.array([0., 0., 0., 0.5, 0.75, 0.75, 1., 1., 1.])
    C = [[ 6.0, 0.0, 6.0],
         [-5.5, 0.2, 5.5],
         [-5.0, 1.3,-5.0],
         [ 4.5, 2.,-4.5],
         [ 4.0, 2.5, 4.0],
         [-3.5, 2.75, 3.5],]

    crv = cad_nurbs([T], C)
    return crv
# ...

# ... custom plot curve
def curve_plot(crv, text="P", color="r" \
               , xmin=None, xmax=None \
               , ymin=None, ymax=None \
               , n_points=100):

    frame1 = plt.gca()
    frame1.axes.get_xaxis().set_visible(False)
    frame1.axes.get_yaxis().set_visible(False)
    frame1.set_frame_on(False)

    if (xmin is not None) and (xmax is not None):
        plt.xlim(xmin, xmax)

    if (ymin is not None) and (ymax is not None):
        plt.ylim(ymin, ymax)

    t = np.linspace(0.0,1.0,n_points)
    C = crv.points
    P = crv(t)
    x = P[:,0]
    y = P[:,1]

    plt.plot(x,y, color)

    i = 0
    for [xc,yc,zc] in C:
        label = "$"+text+"_{%s}$" % (str(i+1))
        plt.text(xc-0.05,yc+0.05, label, fontsize=14, color="k")
        plt.plot(xc,yc, "ok")
        i += 1
    # control polygon
    plt.plot(C[...,0], C[...,1], "-.k")
# ...


# ... create a 3D curve
crv0 = make_curve()
crv1 = make_curve()
# ...

# ... refine the curve by inserting the knot 0.5
crv1.refine(0, [0.5])
# ...

# ... raising the spline degree
crv1.elevate(0,1)
# ...

# ... define the matplotlib grid
xmin = -6.5 ; xmax = 7.
ymin = -0.5 ; ymax = 3.25
# ...

curve_plot(crv0, text="P", color="b")
plt.savefig("curve_ex1_original.png")
plt.clf()

curve_plot(crv1, text="Q", color="r")
plt.savefig("curve_ex1_refined.png")
plt.clf()

The original curve is

_images/curve_ex1_original.png

After inserting a knot and elevating the spline degree, we get

_images/curve_ex1_refined.png

cad_geometry

In this section, we show some useful functions of the cad_geometry object.

Polar Extrude

Let’s start with a 2D curve as given by the following plot

# ... create a 3D curve
def make_curve():
    T = np.array([0., 0., 0., 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1., 1., 1.])
    C = [[ 0.,-2., 0.],
         [ 1.,-2., 0.],
         [ 1., 0., 0.],
         [ 1., 1., 0.],
         [ 0., 1., 0.],
         [-1., 1., 0.],
         [-1., 0., 0.],
         [-1.,-2., 0.],
         [ 0.,-2., 0.],]

    crv = cad_nurbs([T], C)
    return crv
# ...

# ... create a 3D curve
crv0 = make_curve()
geo = cad_geometry()
geo.append(crv0)
geo.plotMesh(MeshResolution=20)
plt.savefig("polar_extrude_ex1_original.png")
plt.clf()
# ...

The result is the following plot

_images/polar_extrude_ex1_original.png

Now we can apply the polar extrude algorithm to fill in the domain

# ... polar extrude
geo_1 = geo.polarExtrude()
U = np.linspace(0., 1., 9)[1:-1]
V = np.linspace(0., 1., 17)[1:-1]
V = [v for v in V if v not in crv0.knots[0]]
geo_1.refine(list_t=[U, V])
geo_1.plotMesh()
plt.savefig("polar_extrude_ex1_geo1.png")
plt.clf()
# ...

The result is the following plot, which is singular at the center.

_images/polar_extrude_ex1_geo1.png

The user can also specify a center, or give a scaling factor, to scale the boundary with respect to the center. The following script shows how to do this

# ... polar extrude
geo_2 = geo.polarExtrude(t=0.25)
U = np.linspace(0., 1., 9)[1:-1]
V = np.linspace(0., 1., 17)[1:-1]
V = [v for v in V if v not in crv0.knots[0]]
geo_2.refine(list_t=[U, V])
geo_2.plotMesh()
plt.savefig("polar_extrude_ex1_geo2.png")
plt.clf()
# ...

The result is the following plot

_images/polar_extrude_ex1_geo2.png

Now, what we want is to fill in the hole and avoid the singular map. Let’s first get some information about the boundary, specially, we need to know what is the face id of the internal face. This can be done by calling

# ... polar extrude into 5 patchs
geo_3 = geo.polarExtrude(t=0.25)
#     print boundaries info in order to know what is the internal boundary
nrb = geo_3[0]
nrb.plotBoundariesInfo()
plt.savefig("polar_extrude_ex1_geo3_info.png")
plt.clf()

The result is the following plot

_images/polar_extrude_ex1_geo3_info.png

Therefor, we use the face id, by calling the to5patchs functions as the following

#     now that we know that the internal face is 1
#     we can use it to split the geometry into 5 patchs
geo_3 = geo_3.to5patchs(1)

U = np.linspace(0., 1., 9)[1:-1]
geo_3.refine(list_t=[U, U])
geo_3.plotMesh()
plt.savefig("polar_extrude_ex1_geo3.png")
plt.clf()
# ...

The result is the following plot

_images/polar_extrude_ex1_geo3.png

Multi-Patchs geometries

In many applications, we need to use a multi-patch description. This means that the physical domain will be splitted into sub-domains, each of them is an image of a square in 2D (line in 1D, cube in 3D). Therefore, we need to provide the connectivity of these macro-elements. For this purpose, we have added the merge function that does this (almost) automatically. Next we give two examples of constructing a 4 patchs description of the square and a 5 patchs description of the circle.

# ... Import the square domain -A-
geo_1 = square(n=[nx,ny],p=[px,py])

# ... Import the square domain -B-
geo_2 = square(n=[nx,ny],p=[px,py])
geo_2[0].translate(1.0, axis=0)

# ... Import the square domain -C-
geo_3 = square(n=[nx,ny],p=[px,py])
geo_3[0].translate(1.0, axis=1)

# ... Import the square domain -D-
geo_4 = square(n=[nx,ny],p=[px,py])
geo_4[0].translate(1.0, axis=0)
geo_4[0].translate(1.0, axis=1)

# ... merging geometries
geo_12 = geo_1.merge(geo_2)
geo_34 = geo_3.merge(geo_4)
geo = geo_12.merge(geo_34)

# ... scaling to the unit square
for i in range(0, geo.npatchs):
    geo[i].scale(0.5, axis=0)
    geo[i].scale(0.5, axis=1)

which gives the following subdivision

_images/square_4mp.png

An interesting example is the construction of a circle domain using 5 patchs, which can be done as the following

# ... Import the quart_circle domain -A-
geo_1 = quart_circle(n=[nx,ny],p=[px,py])
geo_1[0].transpose()

# ... Import the quart_circle domain -B-
geo_2 = quart_circle(n=[nx,ny],p=[px,py])
geo_2[0].rotate(0.5*np.pi)
geo_2[0].transpose()

# ... Import the quart_circle domain -C-
geo_3 = quart_circle(n=[nx,ny],p=[px,py])
geo_3[0].rotate(np.pi)
geo_3[0].reverse(0)

# ... Import the quart_circle domain -D-
geo_4 = quart_circle(n=[nx,ny],p=[px,py])
geo_4[0].rotate(1.5*np.pi)
geo_4[0].reverse(0)

# ... Import the circle domain -E-
geo_5 = circle(radius=0.5,n=[nx,ny],p=[px,py])
geo_5[0].rotate(0.25*np.pi)
geo_5[0].rotate(0.5*np.pi)

# ... merging geometries
geo_12   = geo_1.merge(geo_2)
geo_34   = geo_3.merge(geo_4)
geo_1234 = geo_12.merge(geo_34)
geo      = geo_1234.merge(geo_5)

which gives the following subdivision

_images/circle_5mp.png

Note

Sometimes you need to reverse the parametrization of the domain in order to merge it with another one. In the previous example, we see that we had to reverse the parametrization for the quart circles C and D. Otherwise, pigasus will not reconize the duplicated faces. In this case, you will get a message informing you that pigasus found 2 faces that can match if you reverse the parametrization of one of them.

The resulting mesh for the previous geometry is

_images/circle_5mp_mesh.png

We can check the Jacobian of the corresponding mappings by typing

import matplotlib.pyplot as plt
geo.plotJacobians(MeshResolution=50)
plt.show()

which gives the following plot

_images/circle_5mp_jacobian.png

About the XML format

Many formats exist to describe B-splines and NURBS. We prefered to developp our own format so that we can include some additional data needed for the Finite Element solver.

In some cases, the user may want to add specific informations about one single B-spline/NURBS (patch). This can be achieved using XML attributs. For instance, any cad_object object offers the set_attribut function. Here are some examples:

C = [[0, 1], [1, 1], [1, 0]] # 3x2 grid of 2D control points
w = [1, np.sqrt(2)/2, 1]     # rational weigths
U = [0,0,0, 1,1,1]           # knot vector
crv = cad_nurbs([U], C, weights=w)
crv.set_attribut("name", "arc-circle")
crv.set_attribut("color", "#7A3030")
crv.set_attribut("type", "nurbs")

For the moment only three attributs exist:

  • name : sets a name for the patch,
  • color : sets a color for the patch,
  • type : the type of the patch. Possible values are spline, nurbs and gradiant. When not specified, this means that the patch is either a spline or a nurbs.

IO drivers

Namelist format

ASCII

HDF5

Todo

Implement the hdf5 IO driver

Bezier-Description format

ASCII

HDF5

Todo

Implement the hdf5 IO driver

Hermite-Bezier format

ASCII

HDF5

Todo

Implement the hdf5 IO driver

GeoPDEs format

ASCII