Uduh at 11:32, 1 September 2020

2020-09-01T11:32:02Z

Uduh: Created page with "Go back to Examples. Medusa supports filling domains bounded by non-uniform rational basis splines ([https://en.wikipedia.org/wiki/Non-uniform_rational_B-..."

2020-09-01T11:26:46Z

Created page with "Go back to Examples. Medusa supports filling domains bounded by non-uniform rational basis splines ([https://en.wikipedia.org/wiki/Non-uniform_rational_B-..."

New page

Go back to [[Medusa#Examples|Examples]].

Medusa supports filling domains bounded by non-uniform rational basis splines ([https://en.wikipedia.org/wiki/Non-uniform_rational_B-spline, NURBS]) curves and NURBS surfaces obtained as the tensor product of two NURBS curves. Note that [https://en.wikipedia.org/wiki/B%C3%A9zier_curve Bézier curves] and [https://en.wikipedia.org/wiki/B-spline B-splines] are special cases of NURBS curves.

See [[Positioning of computational nodes]] for details on node positioning algorithms.

== Poisson equation with Dirichlet boundary conditions in 2D ==
With Medusa, we can also solve partial differential equations on NURBS domains. Consider the solution of a simple 2D Poisson equation with Dirichlet boundary conditions:
<math>
\begin{align*}
\Delta u &= 20 &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
</math>
where $\Omega$ is the interior of a NURBS curve defined piecewise. In this case, we will use a duck shaped curve composed of $8$ NURBS curve patches.

In Medusa, we represent NURBS patches as instances of <code>NURBSPatch<vec_t, param_vec_t></code> class, where <code>vec_t</code> is the vector type of the ambient space (in this case <code>Vec2d</code>) and <code>param_vec_t</code> is the vector type of the parametric space (<code>Vec1d</code> for curves and <code>Vec2d</code> for surfaces). To construct an instance, its control points, weights, knot vector, and order $p$ must be given. In this case, we have a simple knot vector <code>{0, 0, 0, 0, 1, 1, 1, 1}</code>, all weights equal to $1$ and $p = 3$ (this means that our curves are Bézier curves). Control points will be calculated from <code>pts</code>. We will push all patches into one <code>Range</code>
<syntaxhighlight lang="c++" line>
// Points to be used in NURBS creation.
Range<Vec2d> pts{{210, -132.5}, {205, -115}, {125, -35}, {85, -61}, {85, -61}, {85, -61},
{80, -65}, {75, -58}, {75, -58}, {65, 45}, {25, 25}, {-15, -16}, {-15, -16},
{-15, -16}, {-35, -28}, {-40, -32.375}, {-40, -32.375}, {-43, -35}, {-27, -40},
{-5, -40}, {-5, -40}, {50, -38}, {35, -65}, {15, -75}, {15, -75}, {-15, -95},
{-20, -140}, {45, -146}, {45, -146}, {95, -150}, {215, -150}, {210, -132.5}};

// Calculate NURBS patches' control points, weights, knot vector and order.
int p = 3;
Range<double> weights{1, 1, 1, 1};
Range<double> knot_vector{0, 0, 0, 0, 1, 1, 1, 1};
Range<NURBSPatch<Vec2d, Vec1d>> patches;

for (int i = 0; i < 8; i++) {
Range<Vec2d> control_points;
for (int j = 0; j < 4; j++) {
control_points.push_back(pts[4 * i + j]);
}

NURBSPatch<Vec2d, Vec1d> patch(control_points, weights, {knot_vector}, {p});
patches.push_back(patch);
}
</syntaxhighlight>

Now, we can construct a <code>NURBSShape<vec_t, param_vec_t></code> representing our NURBS curve from a <code>Range</code> of patches and fill it according to the spacing $h$. Spacing can also be given in the form of a spacing function, see [[Parametric domains]].
<syntaxhighlight lang="c++" line>
// Create NURBS shape from a range of patches.
NURBSShape<Vec2d, Vec1d> shape(patches);

// Fill the domain.
double h = 0.5;
DomainDiscretization<Vec2d> domain = shape.discretizeBoundaryWithStep(h);
</syntaxhighlight>

After that, we can fill the interior of the domain in the usual way.
<syntaxhighlight lang="c++" line>
GeneralFill<Vec2d> gf;
gf(domain, h);
</syntaxhighlight>

Finally, we translate the partial differential equations of our problem into code and solve the problem, as in other examples.
<syntaxhighlight lang="c++" line>
// Construct the approximation engine.
int m = 2; // basis order
Monomials<Vec2d> mon(m);
RBFFD<Polyharmonic<double, 3>, Vec2d> approx({}, mon);

// Find support for the nodes.
int N = domain.size();
domain.findSupport(FindClosest(mon.size()));

// Compute the shapes (we only need the Laplacian).
auto storage = domain.computeShapes<sh::lap>(approx);

Eigen::SparseMatrix<double, Eigen::RowMajor> M(N, N);
Eigen::VectorXd rhs(N);
rhs.setZero();
M.reserve(storage.supportSizes());

// Construct implicit operators over our storage.
auto op = storage.implicitOperators(M, rhs);

for (int i : domain.interior()) {
// set the case for nodes in the domain
op.lap(i) = 20.0;
}
for (int i : domain.boundary()) {
// enforce the boundary conditions
op.value(i) = 0.0;
}

Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
solver.compute(M);
ScalarFieldd u = solver.solve(rhs);
</syntaxhighlight>

Here is the plot of $u(x, y)$.

TODO!!

The whole example can be found as [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation poisson_NURBS_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation poisson_NURBS_2D_2D.m].

@@ Line 15: / Line 15: @@
 where $\Omega$ is the interior of a NURBS curve defined piecewise. In this case, we will use a duck shaped curve composed of $8$ NURBS curve patches.
-In Medusa, we represent NURBS patches as instances of <code>NURBSPatch<vec_t, param_vec_t></code> class, where <code>vec_t</code> is the vector type of the ambient space (in this case <code>Vec2d</code>) and <code>param_vec_t</code> is the vector type of the parametric space (<code>Vec1d</code> for curves and <code>Vec2d</code> for surfaces). To construct an instance, its control points, weights, knot vector, and order $p$ must be given. In this case, we have a simple knot vector <code>{0, 0, 0, 0, 1, 1, 1, 1}</code>, all weights equal to $1$ and $p = 3$ (this means that our curves are Bézier curves). Control points will be calculated from <code>pts</code>. We will push all patches into one <code>Range</code>
+In Medusa, we represent NURBS patches as instances of <code>NURBSPatch<vec_t, param_vec_t></code> class, where <code>vec_t</code> is the vector type of the ambient space (in this case <code>Vec2d</code>) and <code>param_vec_t</code> is the vector type of the parametric space (<code>Vec1d</code> for curves and <code>Vec2d</code> for surfaces). To construct an instance, its control points, weights, knot vector, and order $p$ must be given. In this case, we have a simple knot vector <code>{0, 0, 0, 0, 1, 1, 1, 1}</code>, all weights equal to $1$ and $p = 3$ (this means that our curves are Bézier curves). Control points will be calculated from <code>pts</code>. We will push all patches into one <code>Range</code>.
 <syntaxhighlight lang="c++" line>
 // Points to be used in NURBS creation.
@@ Line 95: / Line 95: @@
 Here is the plot of $u(x, y)$.
-TODO!!
+[[File:duck_pde.png|800px]]
 The whole example can be found as [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation poisson_NURBS_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation poisson_NURBS_2D_2D.m].

NURBS domains - Revision history

Uduh at 11:32, 1 September 2020

Uduh: Created page with "Go back to Examples. Medusa supports filling domains bounded by non-uniform rational basis splines ([https://en.wikipedia.org/wiki/Non-uniform_rational_B-..."