Meshing surface of a sphere with a subdomain
I am trying to build a triangle mesh of the surface of a sphere which also includes a subdomain defined by a 'polygon'. Here is a successful example (subdomain defined by the red dots):
Note that the triangles crossing lon = 180 are not shown.
The method I used to do this fails for finer meshes, because it only 'knows' about the individual points and does not 'know' about the polygon boundary, and producing elements which do not conform to the boundary:
It appears that there are many libraries for meshing in the plane, but I have struggled to find any libraries capable of meshing implicit spherical surfaces with explicit subdomains. Please could anybody suggest a library with this capability? Preferably with a Python API.
For completeness, here is my code:
``` import os
import matplotlib.pyplot as plt import meshio import numpy as np import pygmsh
def rlonlat_to_xyz(r, lon, lat): ''' Converts from radius, longitude and latitude to Cartesian coordinates. '''
theta = (np.pi/2.0) - lat
x = r*np.sin(theta)*np.cos(lon)
y = r*np.sin(theta)*np.sin(lon)
z = r*np.cos(theta)
return x, y, z
def xyz_to_rlonlat(x, y, z): ''' Converts from Cartesian coordinates to radius, longitude and latitude. '''
r = np.sqrt(x**2.0 + y**2.0 + z**2.0)
theta = np.arccos(z/r)
lat = (np.pi/2.0) - theta
lon = np.arctan2(y, x)
return r, lon, lat
def get_real_points(mesh): ''' The gmsh library adds extra points to help in the meshing process. 506
These are returned by generate_mesh(). It is not clear from the pygmsh documentation how to remove these points. Here we remove them by selecting only the points belonging to triangles. It is also necessary to update the triangulation indices. '''
# Get a list of the indices of points used in the triangulation.
tri = mesh.get_cells_type('triangle')
i_mesh = np.sort(np.unique(tri.flatten()))
# Get the points used in the triangulation.
pts = mesh.points[i_mesh, :]
# Define a mapping such that tri_new[i, j] = i_mapping[tri[i, j].
# Note that i_mapping[k] == -1 if k does not belong to i_mesh.
i_max = np.max(i_mesh)
i_mapping = np.zeros((i_max + 1), dtype = np.int) - 1
j = 0
for i in range(i_max + 1):
if i in i_mesh:
i_mapping[i] = j
j = j + 1
# Apply the mapping to create a new triangulation.
n_tri = tri.shape[0]
tri_new = np.zeros(tri.shape, dtype = np.int) - 1
for i in range(n_tri):
for j in range(3):
tri_new[i, j] = i_mapping[tri[i, j]]
return pts, tri_new
def main():
# Set mesh size.
mesh_size = 0.2
#mesh_size = 0.1
# Define polygon longitude and latitude.
lon_pts, lat_pts = np.array([
[130.04592083, <PHONE_NUMBER>, <PHONE_NUMBER>, 152.00842172, <PHONE_NUMBER>,
<PHONE_NUMBER>, <PHONE_NUMBER>, 110.82292747, <PHONE_NUMBER>, <PHONE_NUMBER>,
<PHONE_NUMBER>, <PHONE_NUMBER>, <PHONE_NUMBER>, <PHONE_NUMBER>, <PHONE_NUMBER>,
<PHONE_NUMBER>, <PHONE_NUMBER>, 58. 988