# Mesh adaptation
**URL:** https://pyfr.discourse.group/t/mesh-adaptation/1206
**Category:** Just Starting
**Created:** [15 January 2025 13:58 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206 "2025-01-15T13:58:21Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@mschouler](https://pyfr.discourse.group/u/mschouler)
#### Post date: [15 January 2025 13:58 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/1 "2025-01-15T13:58:21Z")
Hello,
I am currently trying to couple pyfr with the mesh adaptation software [MAdLib](https://sites.uclouvain.be/madlib/) which performs anisotropic adaptation with tri. / tet. of arbitrary order.
So far I have been able to adapt the mesh externally with the solution produced at a given time step but then, the solution needs to be interpolated on the adapted mesh brefore restarting the computation. I was hence wondering if pyfr provides methods to do such interpolation.
Please note that as advised by Peter, I’m currently working off of the HEAD of the develop branch to stay up to date with the latest format developments and I am using the 2D Euler vortex as an example.
Regards,
Marc
---
### Author:  [@fdw](https://pyfr.discourse.group/u/fdw)
#### Post date: [15 January 2025 14:22 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/2 "2025-01-15T14:22:51Z")
Yes; the develop version of PyFR contains support for resampling. See `pyfr resample`. Performance is slightly below where I would like it to be although this is being worked on. Similarly, the algorithm itself isn’t perfect but is reasonable.
Regards, Freddie.
---
### Author:  [@mschouler](https://pyfr.discourse.group/u/mschouler)
#### Post date: [15 January 2025 15:35 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/3 "2025-01-15T15:35:06Z")
Hi Freddie,
I had to subtract 1 to `q[self.comm.rank]` in [`_splitters`](https://github.com/PyFR/PyFR/blob/2db4c41f5f0670c7da34cfea247d62ef4ff7b2c7/pyfr/mpiutil.py#L361) to avoid an overflow error and it seems to work just fine, thank you very much!
Cheers,
Marc
---
### Author:  [@fdw](https://pyfr.discourse.group/u/fdw)
#### Post date: [15 January 2025 16:52 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/4 "2025-01-15T16:52:17Z")
Can you let me know if the following resolves the issue?
```auto
diff --git a/pyfr/mpiutil.py b/pyfr/mpiutil.py
index ba80a3c0..1ef007e7 100644
--- a/pyfr/mpiutil.py
+++ b/pyfr/mpiutil.py
@@ -319,7 +319,7 @@ class SparseScatterer(AlltoallMixin):
class Sorter(AlltoallMixin):
- def __init__ (self, comm, keys, dtype=int):
+ def __init__ (self, comm, keys, dtype=np.uint64):
self.comm = comm
# Locally sort our outbound keys
@@ -354,7 +354,7 @@ class Sorter(AlltoallMixin):
W = math.ceil(np.log2(kmax - kmin + 1))
e, rt = 0, 0
- q = np.empty(self.comm.size, dtype=int)
+ q = np.empty(self.comm.size, dtype=np.uint64)
for i in range(W - 1, -1, -1):
# Compute and gather the probes
diff --git a/pyfr/nputil.py b/pyfr/nputil.py
index 14af796d..f25ff6e0 100644
--- a/pyfr/nputil.py
+++ b/pyfr/nputil.py
@@ -56,7 +56,7 @@ def clean(origfn=None, tol=1e-10, ckwarg='clean'):
return cleanfn(origfn) if origfn else cleanfn
-def morton_encode(ipts, imax, dtype=int):
+def morton_encode(ipts, imax, dtype=np.uint64):
# Allocate the codes
codes = np.zeros(len(ipts), dtype=dtype)
diff --git a/pyfr/resamplers.py b/pyfr/resamplers.py
index 6b25619f..b10c0455 100644
--- a/pyfr/resamplers.py
+++ b/pyfr/resamplers.py
@@ -105,7 +105,7 @@ class BaseCloudResampler(AlltoallMixin):
# Normalise our points and compute their Morton codes
fact = 2**21 if len(pmin) == 3 else 2**32
- ipts = (fact*((pts - pmin) / (pmax - pmin).max())).astype(int)
+ ipts = (fact*((pts - pmin) / (pmax - pmin).max())).astype(np.uint64)
keys = morton_encode(ipts, [fact]*len(pmin))
# Use these codes to sort our points
```
Regards, Freddie.
---
### Author:  [@mschouler](https://pyfr.discourse.group/u/mschouler)
#### Post date: [16 January 2025 09:09 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/5 "2025-01-16T09:09:12Z")
Hi Freddie,
Unfortunately no, it does not.
It first breaks in [`nputil.py`](https://github.com/PyFR/PyFR/blob/2db4c41f5f0670c7da34cfea247d62ef4ff7b2c7/pyfr/nputil.py#L81):
```py
for p, pops in zip((ipt >> ishift).T, ops):
```
as `ishift` must also be typed as `np.uint64` and then an overflow error occurs in [`mpiutil.py`](https://github.com/PyFR/PyFR/blob/2db4c41f5f0670c7da34cfea247d62ef4ff7b2c7/pyfr/mpiutil.py#L100):
```py
v = np.array([v], dtype=dtype)
```
I have tried to propagate the type change further naïvely but without success.
I hope this helps,
Marc
---
### Author:  [@fdw](https://pyfr.discourse.group/u/fdw)
#### Post date: [16 January 2025 22:55 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/6 "2025-01-16T22:55:10Z")
Thanks for working through it. Here is an updated patch.
```auto
diff --git a/pyfr/mpiutil.py b/pyfr/mpiutil.py
index 64e56843..305d5fd3 100644
--- a/pyfr/mpiutil.py
+++ b/pyfr/mpiutil.py
@@ -320,12 +320,17 @@ class SparseScatterer(AlltoallMixin):
class Sorter(AlltoallMixin):
- def __init__ (self, comm, keys, dtype=int):
+ typemap = {
+ np.int8: np.uint8, np.int16: np.uint16,
+ np.int32: np.uint32, np.int64: np.uint64
+ }
+
+ def __init__ (self, comm, keys):
self.comm = comm
# Locally sort our outbound keys
self.sidx = np.argsort(keys)
- skeys = keys[self.sidx].astype(dtype)
+ skeys = keys[self.sidx]
# Determine the total size of the array
size = scal_coll(comm.Allreduce, len(keys))
@@ -346,16 +351,30 @@ class Sorter(AlltoallMixin):
self.ridx = np.argsort(rkeys)
self.keys = rkeys[self.ridx]
+ def _transform_keys(self, skeys):
+ dtype = skeys.dtype
+
+ if np.issubdtype(dtype, np.unsignedinteger):
+ return skeys
+ elif np.issubdtype(dtype, np.signedinteger):
+ tdtype = self.typemap[dtype]
+ return skeys.view(tdtype) ^ tdtype(np.iinfo(dtype).max + 1)
+ else:
+ raise ValueError('Unsupported dtype')
+
def _splitters(self, skeys, r):
+ # Transform the keys so they're unsigned integers
+ skeys = self._transform_keys(skeys)
+
# Determine the minimum and maximum values in the array
- kmin = scal_coll(self.comm.Allreduce, skeys[0], op=mpi.MIN)
- kmax = scal_coll(self.comm.Allreduce, skeys[-1], op=mpi.MAX)
+ kmin = self.comm.allreduce(int(skeys[0]), op=mpi.MIN)
+ kmax = self.comm.allreduce(int(skeys[-1]), op=mpi.MAX)
# Compute the number of bits in the key space
- W = math.ceil(np.log2(kmax - kmin + 1))
+ W = math.ceil(math.log2(kmax - kmin + 1))
e, rt = 0, 0
- q = np.empty(self.comm.size, dtype=int)
+ q = np.empty(self.comm.size, dtype=skeys.dtype)
for i in range(W - 1, -1, -1):
# Compute and gather the probes
@@ -384,7 +403,7 @@ class Sorter(AlltoallMixin):
q[self.comm.rank] = r - rt
self.comm.Allgather(mpi.IN_PLACE, q)
- return lbnd + np.maximum(0, np.minimum(ld, q - gd))
+ return lbnd + np.maximum(0, np.minimum(ld, q.astype(int) - gd))
def __call__ (self, svals):
# Locally sort our data
diff --git a/pyfr/nputil.py b/pyfr/nputil.py
index 14af796d..d26140b2 100644
--- a/pyfr/nputil.py
+++ b/pyfr/nputil.py
@@ -56,14 +56,14 @@ def clean(origfn=None, tol=1e-10, ckwarg='clean'):
return cleanfn(origfn) if origfn else cleanfn
-def morton_encode(ipts, imax, dtype=int):
+def morton_encode(ipts, imax, dtype=np.uint64):
# Allocate the codes
codes = np.zeros(len(ipts), dtype=dtype)
# Determine how many bits to use for each input dimension
ndims = ipts.shape[1]
obits = 8*codes.dtype.itemsize
- ibits = obits // ndims
+ ibits = dtype(obits // ndims)
ishift = np.array([max(int(p).bit_length() - ibits, 0) for p in imax])
# Compute the masks and shifts
diff --git a/pyfr/resamplers.py b/pyfr/resamplers.py
index 456a20d6..b7591373 100644
--- a/pyfr/resamplers.py
+++ b/pyfr/resamplers.py
@@ -105,7 +105,7 @@ class BaseCloudResampler(AlltoallMixin):
# Normalise our points and compute their Morton codes
fact = 2**21 if len(pmin) == 3 else 2**32
- ipts = (fact*((pts - pmin) / (pmax - pmin).max())).astype(int)
+ ipts = (fact*((pts - pmin) / (pmax - pmin).max())).astype(np.uint64)
keys = morton_encode(ipts, [fact]*len(pmin))
# Use these codes to sort our points
@@ -130,7 +130,7 @@ class BaseCloudResampler(AlltoallMixin):
bl = int(keys.max()).bit_length()
# With this break our points up into their top-level orthants
- sidx = np.searchsorted(keys, np.arange(2**nb) << (bl - nb))
+ sidx = np.searchsorted(keys, np.arange(1, 2**nb) << (bl - nb))
# Allocate global bounding box array
bboxes = np.full((comm.size, 2**nb, 2, self.ndims), np.nan)
```
Regards, Freddie.
---
### Author:  [@mschouler](https://pyfr.discourse.group/u/mschouler)
#### Post date: [17 January 2025 08:54 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/7 "2025-01-17T08:54:39Z")
Hi Freddie,
The patch fixes the bug, thank you very much.
Regards,
Marc
---
### Author:  [@fdw](https://pyfr.discourse.group/u/fdw)
#### Post date: [17 January 2025 15:29 UTC](https://pyfr.discourse.group/t/mesh-adaptation/1206/8 "2025-01-17T15:29:51Z")
Let me know if you run into any other problems with the routine as it is still quite fresh.
Regards, Freddie.