#!/usr/bin/env python
# coding: utf-8

# # Création du maillage (IR cylindrique 3D)

# In[1]:


import gmsh
import sys
import math
import numpy

gmsh.initialize(sys.argv)
gmsh.clear()
gmsh.model.add("Test_3D")
gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.5);
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 1);


# ### Définition du sol (troué)

# In[2]:


box = gmsh.model.occ.add_box(0, 0, 0, 1, 1, 1)
cyl = gmsh.model.occ.addCylinder(0.5, 0.5, 0.25 ,0, 0, 0.75, 0.2)
sol = gmsh.model.occ.cut([(3,box)],[(3,cyl)])


# In[3]:


gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
assert(volumes == sol[0])
gmsh.model.addPhysicalGroup(3, [volumes[0][1]], 102)


# ### Définition IR

# In[4]:


cyl2 = gmsh.model.occ.addCylinder(0.5, 0.5, 0.25 ,0, 0, 0.75, 0.2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(3, [cyl2], 101)


# In[5]:


gmsh.option.setNumber("Mesh.MeshSizeMax", 0.08)

gmsh.model.mesh.generate(3)

gmsh.write("Test_3D_mesh.vtk")
gmsh.write("Test_3D_mesh.msh")

gmsh.finalize()


# In[6]:


import pyvista
pyvista.set_plot_theme("document")

p = pyvista.Plotter(window_size=(800, 800))
p.add_mesh(mesh=pyvista.read("Test_3D_mesh.msh"),
           stitle="Legend",
           show_scalar_bar=True,
           show_edges=True)
p.add_text("Sample")
p.view_isometric()
p.show()


# # Calcul sous getfem

# ### Initialisation, définition des sous-domaines

# In[7]:


import getfem as gf
import numpy as np
import math


# In[9]:


m = gf.Mesh('import', 'gmsh','Test_3D_mesh.msh')
m.region(101)


# In[10]:


m.region(102)