Commit 38bd0c4b authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

merging

parents 595abfe2 9f404435
Pipeline #5474 failed with stage
in 19 seconds
......@@ -20,7 +20,7 @@
# see <http://www.gnu.org/licenses/>.
mbtests :
image : immc/migflow-valid:v0.6
image : immc/migflow-valid:v0.7
script:
- mkdir build
- cd build/
......
from distutils.sysconfig import get_python_inc
import platform
import os
import ycm_core
# These are the compilation flags that will be used in case there's no
# compilation database set (by default, one is not set).
# CHANGE THIS LIST OF FLAGS. YES, THIS IS THE DROID YOU HAVE BEEN LOOKING FOR.
flags = [
'-Wno-unused-function',
'-std=gnu99',
'-Iscontact',
'-Ihxt',
'-Isrc',
'-DDIMENSION=2'
]
# Clang automatically sets the '-std=' flag to 'c++14' for MSVC 2015 or later,
# which is required for compiling the standard library, and to 'c++11' for older
# versions.
#if platform.system() != 'Windows':
# flags.append( '-std=c++11' )
# Set this to the absolute path to the folder (NOT the file!) containing the
# compile_commands.json file to use that instead of 'flags'. See here for
# more details: http://clang.llvm.org/docs/JSONCompilationDatabase.html
#
# You can get CMake to generate this file for you by adding:
# set( CMAKE_EXPORT_COMPILE_COMMANDS 1 )
# to your CMakeLists.txt file.
#
# Most projects will NOT need to set this to anything; you can just change the
# 'flags' list of compilation flags. Notice that YCM itself uses that approach.
compilation_database_folder = ''
if os.path.exists( compilation_database_folder ):
database = ycm_core.CompilationDatabase( compilation_database_folder )
else:
database = None
SOURCE_EXTENSIONS = [ '.cpp', '.cxx', '.cc', '.c', '.m', '.mm' ]
def DirectoryOfThisScript():
return os.path.dirname( os.path.abspath( __file__ ) )
def IsHeaderFile( filename ):
extension = os.path.splitext( filename )[ 1 ]
return extension in [ '.h', '.hxx', '.hpp', '.hh' ]
def FindCorrespondingSourceFile( filename ):
if IsHeaderFile( filename ):
basename = os.path.splitext( filename )[ 0 ]
for extension in SOURCE_EXTENSIONS:
replacement_file = basename + extension
if os.path.exists( replacement_file ):
return replacement_file
return filename
def FlagsForFile( filename, **kwargs ):
# If the file is a header, try to find the corresponding source file and
# retrieve its flags from the compilation database if using one. This is
# necessary since compilation databases don't have entries for header files.
# In addition, use this source file as the translation unit. This makes it
# possible to jump from a declaration in the header file to its definition in
# the corresponding source file.
filename = FindCorrespondingSourceFile( filename )
if not database:
return {
'flags': flags,
'include_paths_relative_to_dir': DirectoryOfThisScript(),
'override_filename': filename
}
compilation_info = database.GetCompilationInfoForFile( filename )
if not compilation_info.compiler_flags_:
return None
# Bear in mind that compilation_info.compiler_flags_ does NOT return a
# python list, but a "list-like" StringVec object.
final_flags = list( compilation_info.compiler_flags_ )
# NOTE: This is just for YouCompleteMe; it's highly likely that your project
# does NOT need to remove the stdlib flag. DO NOT USE THIS IN YOUR
# ycm_extra_conf IF YOU'RE NOT 100% SURE YOU NEED IT.
try:
final_flags.remove( '-stdlib=libc++' )
except ValueError:
pass
return {
'flags': final_flags,
'include_paths_relative_to_dir': compilation_info.compiler_working_dir_,
'override_filename': filename
}
FROM ubuntu:18.10
RUN apt update
RUN apt install -y git python petsc-dev python3-numpy make gmsh cmake g++
RUN apt install -y git python petsc-dev python3-numpy make gmsh cmake g++ python3-distutils
VOLUME ["/etc/gitlab-runner"]
......
# build image
docker build -t immc/MigFlow-valid:v0.2 .
docker build -t immc/migflow-valid:v0.2 .
# push image to docker-hub
docker login
docker push immc/MigFlow-valid
docker push immc/migflow-valid
......@@ -35,7 +35,7 @@ if(PETSC_FOUND)
list(APPEND SRC hxt_linear_system_petsc.c)
endif(PETSC_FOUND)
include_directories(. ../tools ${EXTERNAL_INCLUDES})
include_directories(. ${CMAKE_SOURCE_DIR}/tools ${EXTERNAL_INCLUDES})
add_library(mbfluid2 SHARED ${SRC})
target_link_libraries(mbfluid2 ${EXTERNAL_LIBRARIES} mbtools2)
target_compile_definitions(mbfluid2 PUBLIC "-DDIMENSION=2")
......
......@@ -226,11 +226,10 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double unbnd = unold;
if (wbnd->type == BND_WALL) unbnd = unold/2;
else if(vid>0) unbnd = (unold+unext)/2;
if (unbnd<0) {
for (int id = 0; id < D; ++id) {
f0[U+id] += unbnd*((vid<0 ? 0 : data[vid+id])-u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unbnd*rho/c;
f0[U+id] += (unbnd*(vid<0?0:data[vid+id])-unold*u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unold*rho/c;
}
}
}
......@@ -355,10 +354,11 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
double upstar = up[i]+dt/mass*(contact[i]+mass*G[i]-vol*dp[i]);
double drag = gammastar*(upstar-u[i]/c);
problem->particle_force[ip*D+i] = -drag+G[i]*mass-vol*dp[i];
f[U+i] = -drag;
f[U+i] = -drag;//-vol*dp[i];
}
*dfdu = gammastar/c;
*dfddp = gammastar*dt/mass*vol;
*dfddp = gammastar*dt/mass*vol;//-vol;
}
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc,const double cold, const double *dcold,const double *u_mesh,const double *u_mesh_old, const double *du_mesh_vec, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
......@@ -396,6 +396,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f0[U+i] =
rho*(u[i]-uold[i])/dt
+ c*dp[i]
+ G[i]*rhoreduced*c
+ problem->volume_drag*u[i]*mu;
f00[(U+i)*n_fields+U+i] = rho/dt + problem->volume_drag*mu;
......@@ -1276,6 +1277,130 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
return i;
}
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(problem->node_volume);
Mesh *mesh = problem->mesh;
problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->node_volume[i] = 0;
}
for (int iel = 0; iel < mesh->n_elements; ++iel) {
double dxidx[D][D];
const int *el = &mesh->elements[iel*N_N];
const double detj = mesh_dxidx(mesh, iel, dxidx);
const double vol = node_volume_from_detj(detj);
for (int i = 0; i < N_N; ++i)
problem->node_volume[el[i]] += vol;
}
}
#include "vector.h"
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
int *touched_el_flag = malloc(sizeof(int)*mesh->n_elements);
double L2 = pow2(0.002);
int *nodes = NULL;
int *stack = NULL;
double *nodes_w = malloc(sizeof(double)*mesh->n_nodes);
int *neighbours = malloc(sizeof(int)*mesh->n_elements*(D+1));
for (int i = 0; i< mesh->n_elements; ++i){
touched_el_flag[i] = 0;
for (int d = 0; d < D+1; ++d)
neighbours[i*(D+1)+d] = -1;
}
for (int i = 0; i < mesh->n_interfaces; ++i) {
int *inter = mesh->interfaces + i*4;
if (inter[2] < 0) continue;
neighbours[inter[0]*(D+1)+inter[1]%(D+1)] = inter[2];
neighbours[inter[2]*(D+1)+inter[3]%(D+1)] = inter[0];
}
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 0;
nodes_w[i] = -1;
}
for (int i = 0; i < problem->n_particles; ++i) {
if(problem->particle_element_id[i] == -1) continue;
double *xp = problem->particle_position+i*D;
vectorClear(nodes);
vectorClear(stack);
int p = 0;
double wtot = 0;
*vectorPush(&stack) = problem->particle_element_id[i];
touched_el_flag[problem->particle_element_id[i]] = 1;
while (vectorSize(stack) > p) {
int elid = stack[p];
int close = 0;
for (int j = 0; j < N_N; ++j) {
int jnode = mesh->elements[elid*N_N+j];
int found = nodes_w[jnode] != -1;
/*for (int l = 0; l < vectorSize(nodes); ++l){
if (nodes[l] == jnode) {
found = 1;
break;
}
}
*/
if (found) continue;
double *x = mesh->x + jnode*3;
double d2 = 0;
for (int k = 0; k < D; ++k)
d2 += pow2(x[k]-xp[k]);
*vectorPush(&nodes) = jnode;
if (d2 < L2) {
close = 1;
double w = pow2(sqrt(L2)-sqrt(d2))*problem->node_volume[jnode];
nodes_w[jnode] = w;
wtot += w;
}
else {
nodes_w[jnode] = 0;
}
}
if (close) {
for (int j = 0; j < N_N; ++j) {
int nid = neighbours[elid*(D+1)+j];
if (nid == -1 || touched_el_flag[nid]!=0) continue;
touched_el_flag[nid] = 1;
*vectorPush(&stack) = nid;
}
}
p += 1;
}
if (vectorSize(nodes) >= 3) {
for (size_t j = 0; j < vectorSize(nodes); ++j) {
problem->porosity[nodes[j]] += volume[i]*nodes_w[nodes[j]]/wtot;
nodes_w[nodes[j]] = -1;
}
}
else {
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*D],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j){
problem->porosity[el[j]] += sf[j]*volume[i];
}
}
for (int j = 0; j < vectorSize(stack); ++j) {
touched_el_flag[stack[j]] = 0;
}
}
for (int i = 0; i < mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
problem->porosity[i] = 0.;
}
else{
problem->porosity[i] = 1-problem->porosity[i]/problem->node_volume[i];
}
}
free(neighbours);
vectorFree(nodes);
vectorFree(stack);
free(nodes_w);
}
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)
{
......
......@@ -60,18 +60,17 @@ def _read_array(a) :
def _write_index(idxname,filename,i,t,reset_index) :
if reset_index is None : reset_index = (i==0 or not os.path.isfile(idxname))
if reset_index :
f = open(idxname,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">\n')
f.write(b'<Collection>\n')
else:
f = open(idxname,"rb+")
f.seek(-25,os.SEEK_END)
f.write(b'<DataSet timestep="%.16g" group="" part="0" file="%s"/>\n'%(t,filename.encode()))
f.write(b'</Collection>\n')
f.write(b'</VTKFile>\n')
f.close()
with open(idxname,"wb" if reset_index else "rb+") as f :
if reset_index :
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">\n')
f.write(b'<Collection>\n')
else:
f.seek(-25,os.SEEK_END)
f.write(b'<DataSet timestep="%.16g" group="" part="0" file="%s"/>\n'%(t,filename.encode()))
f.write(b'</Collection>\n')
f.write(b'</VTKFile>\n')
f.close()
def write(basename,i,t,elements,x,data=None,field_data=None,cell_data=None,reset_index=None,vtktype="vtu") :
"""Write VTK UnstructuredMesh (vtu) or PolyData (vtp) files.
......@@ -105,43 +104,43 @@ def write(basename,i,t,elements,x,data=None,field_data=None,cell_data=None,reset
filename = "%s_%05d.%s"%(basename,i,vtktype)
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
vtkdescr = b"UnstructuredGrid" if vtktype=="vtu" else b"PolyData"
f.write(b'<VTKFile type="%s" format="0.1" compressor="vtkZLibDataCompressor">\n' % vtkdescr)
f.write(b'<%s>\n'%vtkdescr)
nel = elements[0].shape[0] if elements is not None else 0
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
_write_array(f,x,b'NumberOfComponents="3"')
f.write(b'</Points>\n')
if elements is not None:
(types,connectivity,offsets) = elements
f.write(b'<Cells>\n')
_write_array(f,connectivity,b'Name="connectivity"')
_write_array(f,offsets,b'Name="offsets"')
_write_array(f,types,b'Name="types"')
f.write(b'</Cells>\n')
if data is not None :
f.write(b'<PointData>\n')
for name,v in data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</PointData>\n')
if cell_data is not None :
f.write(b'<CellData>\n')
for name,v in cell_data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</CellData>\n')
f.write(b'</Piece>\n')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
_write_array(f,v,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,v.shape[0],v.shape[1]))
f.write(b'</FieldData>\n');
f.write(b'</%s>\n' % vtkdescr)
f.write(b'</VTKFile>\n')
f.close()
with open(filename,"wb") as f :
f.write(b'<?xml version="1.0"?>\n')
vtkdescr = b"UnstructuredGrid" if vtktype=="vtu" else b"PolyData"
f.write(b'<VTKFile type="%s" format="0.1" compressor="vtkZLibDataCompressor">\n' % vtkdescr)
f.write(b'<%s>\n'%vtkdescr)
nel = elements[0].shape[0] if elements is not None else 0
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
_write_array(f,x,b'NumberOfComponents="3"')
f.write(b'</Points>\n')
if elements is not None:
(types,connectivity,offsets) = elements
f.write(b'<Cells>\n')
_write_array(f,connectivity,b'Name="connectivity"')
_write_array(f,offsets,b'Name="offsets"')
_write_array(f,types,b'Name="types"')
f.write(b'</Cells>\n')
if data is not None :
f.write(b'<PointData>\n')
for name,v in data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</PointData>\n')
if cell_data is not None :
f.write(b'<CellData>\n')
for name,v in cell_data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</CellData>\n')
f.write(b'</Piece>\n')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
_write_array(f,v,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,v.shape[0],v.shape[1]))
f.write(b'</FieldData>\n');
f.write(b'</%s>\n' % vtkdescr)
f.write(b'</VTKFile>\n')
f.close()
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def read(fname) :
......@@ -158,19 +157,20 @@ def read(fname) :
node_data,cell_data,field_data -- dictionaries whose keys are the name of the node,cell and field data.
"""
root = ET.parse(open(fname,"rb")).getroot()
ftype = root.attrib["type"]
grid = root.find(ftype)
piece = grid.find("Piece")
x = _read_array(piece.find("Points/DataArray"))
cells = piece.find("Cells")
fd = grid.find("FieldData")
cd = piece.find("CellData")
pd = piece.find("PointData")
c = dict((f.attrib["Name"],_read_array(f)) for f in cells) if cells else None
fdata = dict((f.attrib["Name"],_read_array(f)) for f in fd) if fd else None
cdata = dict((f.attrib["Name"],_read_array(f)) for f in cd) if cd else None
pdata = dict((f.attrib["Name"],_read_array(f)) for f in pd) if pd else None
with open(fname,"rb") as f :
root = ET.parse(f).getroot()
ftype = root.attrib["type"]
grid = root.find(ftype)
piece = grid.find("Piece")
x = _read_array(piece.find("Points/DataArray"))
cells = piece.find("Cells")
fd = grid.find("FieldData")
cd = piece.find("CellData")
pd = piece.find("PointData")
c = dict((f.attrib["Name"],_read_array(f)) for f in cells) if cells else None
fdata = dict((f.attrib["Name"],_read_array(f)) for f in fd) if fd else None
cdata = dict((f.attrib["Name"],_read_array(f)) for f in cd) if cd else None
pdata = dict((f.attrib["Name"],_read_array(f)) for f in pd) if pd else None
return x,c,pdata,cdata,fdata
def write_multipart(basename,parts,i,t,reset_index=None):
......@@ -188,15 +188,15 @@ def write_multipart(basename,parts,i,t,reset_index=None):
filename = basename+"_%05d.vtm"%i
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="vtkMultiBlockDataSet" version="1.0" byte_order="LittleEndian">\n');
f.write(b'<vtkMultiBlockDataSet>\n');
for ip,n in enumerate(parts):
f.write(b'<DataSet index="%i" file="%s"/>\n'%(ip,n.encode()));
f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
with open(filename,"wb") as f :
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="vtkMultiBlockDataSet" version="1.0" byte_order="LittleEndian">\n');
f.write(b'<vtkMultiBlockDataSet>\n');
for ip,n in enumerate(parts):
f.write(b'<DataSet index="%i" file="%s"/>\n'%(ip,n.encode()));
f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def string_array_encode(a) :
......
......@@ -302,6 +302,9 @@ class FluidProblem :
self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
nit = self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it))
return nit
def advance_concentration(self,dt):
if self._n_fluids == 2 :
nv = np.linalg.norm(self.solution()[:,:self._dim],axis=1)
nvmax = np.max(nv[self.elements()],axis=1,keepdims=True)
......@@ -312,7 +315,6 @@ class FluidProblem :
print("sub-iterating advection for cfl : %i sub-iterations"%nsub)
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
return nit
def set_particles(self, mass, volume, position, velocity, contact, init=False, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
......
......@@ -57,112 +57,112 @@ class Mesh() :
self.read(filename)
def write(self, filename, format3=None) :
output = open(filename, "w")
if format3 == None :
format3 = self.useFormat3
if format3 :
output.write("$MeshFormat\n3 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Entities\n%i\n" % (len(self.entities)))
for e in self.entities :
output.write("%i %i %i %s\n" % (e.tag, e.dimension, len(e.physicals), " ".join(repr(ip) for ip in e.physicals)))
output.write("$EndEntities\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g 0\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = len(e.vertices) + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
output.write("%i %i %i %i %s %s\n" % (e.tag, e.etype.tag, entity.tag, ntag, " ".join(repr(v[3]) for v in e.vertices), spart))
output.write("$EndElements\n")
else :
output.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = 2 + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
ptag = entity.physicals[0] if entity.physicals else 0
output.write("%i %i %i %i %i %s %s\n" % (e.tag, e.etype.tag, ntag, ptag, entity.tag, spart, " ".join(repr(v[3]) for v in e.vertices)))
output.write("$EndElements\n")
with open(filename, "w") as output :
if format3 == None :
format3 = self.useFormat3
if format3 :
output.write("$MeshFormat\n3 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Entities\n%i\n" % (len(self.entities)))
for e in self.entities :
output.write("%i %i %i %s\n" % (e.tag, e.dimension, len(e.physicals), " ".join(repr(ip) for ip in e.physicals)))
output.write("$EndEntities\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g 0\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = len(e.vertices) + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
output.write("%i %i %i %i %s %s\n" % (e.tag, e.etype.tag, entity.tag, ntag, " ".join(repr(v[3]) for v in e.vertices), spart))
output.write("$EndElements\n")
else :
output.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))