To construct stable and reliable shape functions the support domains need to be non-degenerated <ref name="LeeLocal">Lee CK, Liu X, Fan SC. Local muliquadric approximation for solving boundary value problems. Comput Mech. 2003;30:395-409.</ref>, i.e. the distances between support nodes have to be balanced. Naturally, this condition is fulfilled in regular nodal distributions, but when working with complex geometries, the nodes have to be positioned accordingly. There are different algorithms designed to optimally fill the domain with different shapes <ref name="LohnerGeneral">Löhner R, Oñate E. A general advancing front technique for filling space with arbitrary objects. Int J Numer Meth Eng. 2004;61:1977-91.</ref><ref name="LiuNode">Liu Y, Nie Y, Zhang W, Wang L. Node placement method by bubble simulation and its application. CMES-Comp Model Eng. 2010;55:89.</ref>, see [[Positioning of computational nodes]] for our algorithms. Here, an intrinsic feature of the MLSM is used to take care of that problem. The goal is to minimize the overall support domain degeneration in order to attain stable numerical solution. In other words, a global optimization problem with the overall deformation of the local support domains acting as the cost function is tackled. We seek the global minimum by a local iterative approach. In each iteration, the computational nodes are translated according to the local derivative of the potential
\begin{equation}
\delta \b{p}\left( \b{p} \right)=-\sigma_{k}\sum\limits_{n=1}^{n}{\nabla }V\left( \b{p}-\b{p}_n \right)
\end{equation}
where $V, n, \delta \b{p}, \b{p}_{n}$ and $\sigma_{k}$ stand for the potential, number of support nodes, offset of the node, position of n-th support node and relaxation parameter, respectively. After offsets in all nodes are computed, the nodes are repositioned as
$\b{p}\leftarrow \b{p}+\delta \b{p}\left( \b{p} \right)$.
Presented iterative process procedure begins with positioning of boundary nodes, which is considered as the definition of the domain, and then followed by the positioning of internal nodes.

The <code>BasicRelax</code> engine supports two call types:
* with supplied distribution function <code>relax(domain, func)</code>, where it tries to satisfy the user supplied nodal density function <code>func</code>. This can be achieved only when there is the total number of domain nodes the same as integral of density function over the domain. If there is too much nodes a volatile relax might occur. If there is not enough nodes the relax might become lazy. The best use of this mode is in combination with <code>GeneralFill</code> and/or <code>GeneralSurfaceFill</code>.
* without distribution, where nodes always move towards less populated area. The relax magnitude is simply determined from Annealing factor and distance to the closest node. A simple and stable approach, however, note that this relax always converges towards uniformly distributed nodes.

Example of filling and relaxing a 2D domain can be found in below code snippet and Figures. Note the difference between relax without supplied distribution (<xr id="fig:relax_uniform"/>) and with supplied distribution (<xr id="fig:relax_non-uniform"/>). These examples are extreme, relaxation is usually used to convert an almost desirable node distribution into a desirable node distribution.
More examples can be found in main repository under tests.
<syntaxhighlight lang="c++" line>
// Define domain shape.
double r = 0.25;
BallShape<Vec2d> bs({0.5, 0.5}, r);

// Fill domain disretization with h.
auto h = [](Vec2d p){
return (0.005 + (p[0] - 0.5) * (p[0] - 0.5) / 2 + (p[1] - 0.5) * (p[1] - 0.5) / 2);
};
DomainDiscretization<Vec2d> domain = bs.discretizeBoundaryWithDensity(h);
GeneralFill<Vec2d> gf;
gf(domain, h);

// Relax domain to uniform distribution.
BasicRelax relax;
relax.iterations(1000).initialHeat(1).finalHeat(0).boundaryProjectionThreshold(0.75)
.projectionType(BasicRelax::PROJECT_IN_DIRECTION).numNeighbours(5);
relax(domain);
</syntaxhighlight>

<figure id="fig:relax_uniform">
[[File:relax_uniform.png|thumb|center|800px|<caption> Relaxation of a non-uniform node distribution to a uniform node distribution. </caption>]]
</figure>

<figure id="fig:relax_non-uniform">
[[File:relax_non-uniform.png|thumb|center|800px|<caption> Relaxation of a uniform node distribution to a non-uniform node distribution. </caption>]]
</figure>

=References=
<references/>

Positioning of computational nodes

2020-09-15T12:24:30Z

Uduh:

2020-09-11T11:27:33Z

2020-09-03T08:57:29Z

Uduh: File uploaded with MsUpload

File uploaded with MsUpload

File:3d simple.png

2020-09-03T08:57:29Z

Uduh: File uploaded with MsUpload

File uploaded with MsUpload

Positioning of computational nodes

2020-09-03T08:37:16Z

Uduh:

NURBS domains

2020-09-01T11:32:02Z

Uduh:

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)$.

[[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].