Commit e9b90ced authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

vtk contact output + paraview filter

parent 1211487d
......@@ -959,7 +959,67 @@ void particleProblemWrite(const ParticleProblem *p, const char *filename) {
fclose(output);
}
int particleProblemWriteVtk(ParticleProblem *p, const char *filename) {
static void _writeContactVtkSingle(ParticleProblem *p, FILE *f, int ctype, const char *cname)
{
size_t n = 0;
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
if (p->contacts[i].type == ctype)
n++;
}
fprintf(f, "<DataArray Name=\"%s\" NumberOfTuples=\"%zu\" NumberOfComponents=\"1\" type=\"Float32\" format=\"ascii\">\n",cname, n);
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
if (p->contacts[i].type == ctype) {
fprintf(f, "%.2g ", p->contacts[i].dv);
}
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray Name=\"%s_idx\" NumberOfTuples=\"%zu\" NumberOfComponents=\"2\" type=\"UInt64\" format=\"ascii\">\n",cname,n);
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
if (p->contacts[i].type == ctype) {
fprintf(f, "%zu %zu ", p->contacts[i].o0, p->contacts[i].o1);
}
}
fprintf(f, "\n</DataArray>\n");
}
static int particleProblemWriteContactVtp(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/contacts_%05d.vtp",dirname, id);
FILE *f = fopen(filename, "w");
if (!f){
printf("cannot open file %s\n", filename);
return -1;
}
size_t n_contacts[DIMENSION+2] = {0};
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
n_contacts[p->contacts[i].type]++;
}
fprintf(f, "<VTKFile type=\"PolyData\">\n");
fprintf(f, "<PolyData>\n");
fprintf(f, "<FieldData>\n");
_writeContactVtkSingle(p,f,PARTICLE_PARTICLE,"particle_particle");
_writeContactVtkSingle(p,f,PARTICLE_DISK,"particle_disk");
_writeContactVtkSingle(p,f,PARTICLE_SEGMENT,"particle_segment");
#if DIMENSION == 3
_writeContactVtkSingle(p,f,PARTICLE_TRIANGLE,"particle_triangle");
#endif
fprintf(f, "</FieldData>\n");
fprintf(f, "<Piece NumberOfPoints=\"1\" NumberOfCells=\"0\">\n");
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
fprintf(f, "0 0 0\n");
fprintf(f, "</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "</Piece>\n");
fprintf(f, "</PolyData>\n");
fprintf(f, "</VTKFile>\n");
fclose(f);
free(filename);
return 0;
}
static int particleProblemWriteParticleVtp(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/particles_%05d.vtp",dirname, id);
FILE *f = fopen(filename, "w");
if (!f){
printf("cannot open file %s\n", filename);
......@@ -970,9 +1030,9 @@ int particleProblemWriteVtk(ParticleProblem *p, const char *filename) {
if (p->contacts[i].type == PARTICLE_PARTICLE)
n_pp_contacts ++;
}
fprintf(f, "<VTKFile type=\"UnstructuredGrid\">\n");
fprintf(f, "<UnstructuredGrid>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"%lu\">\n", vectorSize(p->particles), n_pp_contacts);
fprintf(f, "<VTKFile type=\"PolyData\">\n");
fprintf(f, "<PolyData>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"0\">\n", vectorSize(p->particles));
fprintf(f, "<PointData Scalars=\"Radius\" Vectors=\"Velocity\">\n");
fprintf(f, "<DataArray Name=\"Radius\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
......@@ -999,94 +1059,124 @@ int particleProblemWriteVtk(ParticleProblem *p, const char *filename) {
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "<Cells>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n");
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
if (p->contacts[i].type == PARTICLE_PARTICLE) {
fprintf(f, "%i %i ", p->contacts[i].o0, p->contacts[i].o1);
}
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n");
for (size_t i = 0; i < n_pp_contacts; ++i) {
fprintf(f, "%i ", (i+1)*2);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
for (size_t i = 0; i < n_pp_contacts; ++i) {
fprintf(f, "3 ");
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Cells>\n");
fprintf(f, "<CellData Scalars=\"Reaction\">\n");
fprintf(f, "<DataArray Name=\"Reaction\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
if (p->contacts[i].type == PARTICLE_PARTICLE) {
fprintf(f, "%g ", p->contacts[i].dv);
}
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</CellData>\n");
fprintf(f, "</Piece>\n");
fprintf(f, "</UnstructuredGrid>\n");
fprintf(f, "</PolyData>\n");
fprintf(f, "</VTKFile>\n");
fclose(f);
return 0;
}
void particleProblemWriteBoundaryVtk(ParticleProblem *p, const char *filename) {
static int particleProblemWriteBndVtu(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/boundaries_%05d.vtu",dirname, id);
size_t nd = vectorSize(p->disks);
size_t ns = vectorSize(p->segments);
#if DIMENSION == 3
size_t nt = vectorSize(p->triangles);
#else
size_t nt = 0;
#endif
FILE *f = fopen(filename, "w");
f = fopen(filename, "w");
fprintf(f, "<VTKFile type=\"UnstructuredGrid\">\n");
fprintf(f, "<UnstructuredGrid>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfLines=\"%lu\">\n", vectorSize(p->segments)*2, vectorSize(p->segments));
fprintf(f, "<FieldData>\n");
fprintf(f, "<DataArray Name=\"sizes\" NumberOfTuples=\"3\" NumberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n");
fprintf(f, "%i %i %i\n", (int)vectorSize(p->disks), (int)vectorSize(p->segments), 0);
fprintf(f, "</DataArray>\n");
fprintf(f, "</FieldData>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"%lu\">\n", nt*3+ns*2+nd, nt+ns+nd);
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < vectorSize(p->segments); ++i) {
for (size_t i = 0; i < nd; ++i) {
Disk *d = &p->disks[i];
fprintf(f, "%.16g %.16g %.16g ", d->x[0], d->x[1], DIMENSION == 2 ? 0. : d->x[2]);
}
for (size_t i = 0; i < ns; ++i) {
Segment *s = &p->segments[i];
fprintf(f, "%.16g %.16g %.16g ", s->p[0][0], s->p[0][1], DIMENSION == 2 ? 0. : s->p[0][2]);
fprintf(f, "%.16g %.16g %.16g ", s->p[1][0], s->p[1][1], DIMENSION == 2 ? 0. : s->p[1][2]);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "<Lines>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"connectivity\">\n");
for (size_t i = 0; i < vectorSize(p->segments); ++i) {
fprintf(f, "%i %i ", (int)i*2, (int)i*2+1);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"offsets\">\n");
for (size_t i = 0; i < vectorSize(p->segments); ++i) {
fprintf(f, "%i ", (int)(i+1)*2);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Lines>\n");
fprintf(f, "</Piece>\n");
#if DIMENSION == 3
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfPolys=\"%lu\">\n", vectorSize(p->triangles)*3, vectorSize(p->triangles));
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->triangles); ++i) {
for (int i = 0; i < nt; ++i) {
Triangle *t = &p->triangles[i];
fprintf(f, "%.16g %.16g %.16g ", t->p[0][0], t->p[0][1], t->p[0][2]);
fprintf(f, "%.16g %.16g %.16g ", t->p[1][0], t->p[1][1], t->p[1][2]);
fprintf(f, "%.16g %.16g %.16g ", t->p[2][0], t->p[2][1], t->p[2][2]);
}
#endif
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "<Polys>\n");
fprintf(f, "<Cells>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"connectivity\">\n");
for (int i = 0; i < vectorSize(p->triangles)*3; ++i) {
fprintf(f, "%i ", i);
for (size_t i = 0; i < nd; ++i) {
fprintf(f, "%i ", (int)(i));
}
for (size_t i = 0; i < ns; ++i) {
fprintf(f, "%i %i ", (int)(nd+i*2), (int)(nd+i*2+1));
}
#if DIMENSION == 3
for (size_t i = 0; i < nt; ++i) {
fprintf(f, "%i %i %i", (int)(nd+ns*2+i*3+0), (int)(nd+ns*2+i*3+1), (int)(nd+ns*2+i*3+2));
}
#endif
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"offsets\">\n");
for (int i = 0; i < vectorSize(p->triangles); ++i) {
fprintf(f, "%i ", (i+1)*3);
for (size_t i = 0; i < nd; ++i) {
fprintf(f, "%i ", (int)(i+1));
}
for (size_t i = 0; i < ns; ++i) {
fprintf(f, "%i ", (int)nd+(int)(i+1)*2);
}
for (size_t i = 0; i < nt; ++i) {
fprintf(f, "%i ", (int)(nd+ns*2+(i+1)*3));
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Polys>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"types\">\n");
for (size_t i = 0; i < nd; ++i) {
fprintf(f, "%i ", 1);
}
for (size_t i = 0; i < ns; ++i) {
fprintf(f, "%i ", 3);
}
for (size_t i = 0; i < nt; ++i) {
fprintf(f, "%i ", 5);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Cells>\n");
fprintf(f, "</Piece>\n");
#endif
fprintf(f, "</UnstructuredGrid>\n");
fprintf(f, "</VTKFile>\n");
fclose(f);
free(filename);
return 0;
}
int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id) {
if (particleProblemWriteParticleVtp(p,dirname,id) != 0)
return -1;
if (particleProblemWriteBndVtu(p,dirname,id) != 0)
return -1;
if (particleProblemWriteContactVtp(p,dirname,id) != 0)
return -1;
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/particles_problem_%05d.vtm",dirname, id);
FILE *f = fopen(filename, "w");
if (!f)
return -1;
fprintf(f, "<VTKFile type=\"vtkMultiBlockDataSet\" version=\"1.0\" byte_order=\"LittleEndian\">\n");
fprintf(f, "<vtkMultiBlockDataSet>\n");
fprintf(f, "<DataSet index=\"0\" file=\"particles_%05d.vtp\">\n",id);
fprintf(f, "</DataSet>\n");
fprintf(f, "<DataSet index=\"1\" file=\"boundaries_%05d.vtu\">\n",id);
fprintf(f, "</DataSet>\n");
fprintf(f, "<DataSet index=\"2\" file=\"contacts_%05d.vtp\">\n",id);
fprintf(f, "</DataSet>\n");
fprintf(f, "</vtkMultiBlockDataSet>\n");
fprintf(f,"</VTKFile>\n");
fclose(f);
free(filename);
return 0;
}
void particleProblemLoad(ParticleProblem *p, const char *filename)
......
......@@ -27,8 +27,7 @@ double *particleProblemPosition(ParticleProblem *p);
double *particleProblemExternalForces(ParticleProblem *p);
unsigned long int particleProblemNParticle(ParticleProblem *p);
void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax);
int particleProblemWriteVtk(ParticleProblem *p, const char *filename);
void particleProblemWriteBoundaryVtk(ParticleProblem *p, const char *filename);
int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
int *particleProblemSegmentTag(ParticleProblem *p);
......
......@@ -81,14 +81,7 @@ class ParticleProblem :
shutil.move(fname+"_tmp", fname)
def write_vtk(self, odir, i, t) :
outputname = "%s/part-%05i.vtu" % (odir, i)
lib.particleProblemWriteVtk(self._p, (outputname+"_tmp").encode())
shutil.move(outputname + "_tmp", outputname)
def write_boundary_vtk(self, odir, i, t) :
outputname = "%s/boundary-%05i.vtp" % (odir, i)
lib.particleProblemWriteBoundaryVtk(self._p, (outputname+"_tmp").encode())
shutil.move(outputname + "_tmp", outputname)
lib.particleProblemWriteVtk(self._p, odir.encode(), i)
def add_boundary_disk(self, x0, r, tag) :
lib.particleProblemAddBoundaryDisk.restype = c_size_t
......
......@@ -81,13 +81,7 @@ class ParticleProblem :
shutil.move(fname+"_tmp", fname)
def write_vtk(self, odir, i, t) :
outputname = "%s/part-%05i.vtu" % (odir, i)
lib.particleProblemWriteVtk(self._p, (outputname+"_tmp").encode())
shutil.move(outputname + "_tmp", outputname)
def write_boundary_vtk(self, odir, i, t) :
outputname = "%s/boundary-%05i.vtp" % (odir, i)
lib.particleProblemWriteBoundaryVtk(self._p, (outputname+"_tmp").encode())
lib.particleProblemWriteVtk(self._p, odir.encode(), i)
shutil.move(outputname + "_tmp", outputname)
def add_boundary_disk(self, x0, r, tag) :
......
......@@ -672,7 +672,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
if(problem->particle_element_id[i] == -1) continue;
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*DIMENSION],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
const uint32_t *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];
}
......
all : contacts-extract.xml
%.xml : %.py
python python_filter_generator.py $< $@
Name = 'MBExtractContact'
Label = 'Marblesbag extract contacts'
Help = 'Generate tubes to represent contacts from marblesbag scontact output'
NumberOfInput = 1
InputDataType = 'vtkMultiBlockDataSet'
OutputDataType = 'vtkUnstructuredGrid'
ExtraXml = ''
Properties = dict(nf = 5, rf = 2.5e-4)
def RequestData():
import numpy as np
#input
#nf = 5
#rf = 2.5e-4
in_particles,in_bnd,in_contacts = list(inputs[0])
#split bnd cells
n_disks, n_segments, n_triangles = in_bnd.FieldData["sizes"]
if n_disks :
disks = in_bnd.Cells[1:2*n_disks:2]
if n_segments :
segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
if n_triangles :
triangles = in_bnd.Cells[2*n_disks+3*n_segments,2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
vals = []
points = []
#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
vals.append(in_contacts.FieldData["particle_disk"])
contacts = in_contacts.FieldData["particle_disk_idx"]
points_d = np.ndarray((contacts.shape[0],2,3))
points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
points.append(points_d)
#particle-segments contacts
vals.append(in_contacts.FieldData["particle_segment"])
contacts = in_contacts.FieldData["particle_segment_idx"]
points_s = np.ndarray((contacts.shape[0],2,3))
points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
points_s[:,1,:] = np.mean(in_bnd.Points[segments[contacts[:,0],:]],axis=1)
points.append(points_s)
#particle-triangle contacts
if "particle_triangle" in in_contacts.FieldData :
vals.append(in_contacts.FieldData["particle_triangle"])
contacts = in_contacts.FieldData["particle_triangle_idx"]
points_t = np.ndarray((contacts.shape[0],2,3))
points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
points.append(points_t)
#merge everything
points = np.vstack(points)
vals = np.hstack(vals)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(t[:,1,None]>t[:,2,None],np.array([[0,0,1]]),np.array([[0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1
opoints = points[:,None,:,:] \
+n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
+n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
j = (i+1)%nf
pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
<ServerManagerConfiguration>
<ProxyGroup name="filters">
<SourceProxy name="MBExtractContact" class="vtkPythonProgrammableFilter" label="Marblesbag extract contacts">
<Documentation
long_help="Generate tubes to represent contacts from marblesbag scontact output"
short_help="Generate tubes to represent contacts from marblesbag scontact output">
</Documentation>
<InputProperty
name="Input"
command="SetInputConnection">
<ProxyGroupDomain name="groups">
<Group name="sources"/>
<Group name="filters"/>
</ProxyGroupDomain>
<DataTypeDomain name="input_type">
<DataType value="vtkMultiBlockDataSet"/>
</DataTypeDomain>
</InputProperty>
<IntVectorProperty
name="nf"
label="nf"
initial_string="nf"
command="SetParameter"
animateable="1"
default_values="5"
number_of_elements="1">
<Documentation></Documentation>
</IntVectorProperty>
<DoubleVectorProperty
name="rf"
label="rf"
initial_string="rf"
command="SetParameter"
animateable="1"
default_values="0.00025"
number_of_elements="1">
<Documentation></Documentation>
</DoubleVectorProperty>
<!-- Output data type: "vtkUnstructuredGrid" -->
<IntVectorProperty command="SetOutputDataSetType"
default_values="4"
name="OutputDataSetType"
number_of_elements="1"
panel_visibility="never">
<Documentation>The value of this property determines the dataset type
for the output of the programmable filter.</Documentation>
</IntVectorProperty>
<StringVectorProperty
name="Script"
command="SetScript"
number_of_elements="1"
default_values="import numpy as np&#xA;&#xA;#input&#xA;#nf = 5&#xA;#rf = 2.5e-4&#xA;&#xA;in_particles,in_bnd,in_contacts = list(inputs[0])&#xA;&#xA;#split bnd cells&#xA;n_disks, n_segments, n_triangles = in_bnd.FieldData[&quot;sizes&quot;]&#xA;if n_disks :&#xA; disks = in_bnd.Cells[1:2*n_disks:2]&#xA;if n_segments :&#xA; segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]&#xA;if n_triangles :&#xA; triangles = in_bnd.Cells[2*n_disks+3*n_segments,2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]&#xA;&#xA;vals = []&#xA;points = []&#xA;&#xA;&#xA;#particle-particle contacts&#xA;vals.append(in_contacts.FieldData[&quot;particle_particle&quot;])&#xA;contacts = in_contacts.FieldData[&quot;particle_particle_idx&quot;]&#xA;points.append(in_particles.Points[contacts,:])&#xA;&#xA;#particle-disk contacts&#xA;vals.append(in_contacts.FieldData[&quot;particle_disk&quot;])&#xA;contacts = in_contacts.FieldData[&quot;particle_disk_idx&quot;]&#xA;points_d = np.ndarray((contacts.shape[0],2,3))&#xA;points_d[:,0,:] = in_particles.Points[contacts[:,1],:]&#xA;points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]&#xA;points.append(points_d)&#xA;&#xA;#particle-segments contacts&#xA;vals.append(in_contacts.FieldData[&quot;particle_segment&quot;])&#xA;contacts = in_contacts.FieldData[&quot;particle_segment_idx&quot;]&#xA;points_s = np.ndarray((contacts.shape[0],2,3))&#xA;points_s[:,0,:] = in_particles.Points[contacts[:,1],:]&#xA;points_s[:,1,:] = np.mean(in_bnd.Points[segments[contacts[:,0],:]],axis=1)&#xA;points.append(points_s)&#xA;&#xA;#particle-triangle contacts&#xA;if &quot;particle_triangle&quot; in in_contacts.FieldData :&#xA; vals.append(in_contacts.FieldData[&quot;particle_triangle&quot;])&#xA; contacts = in_contacts.FieldData[&quot;particle_triangle_idx&quot;]&#xA; points_t = np.ndarray((contacts.shape[0],2,3))&#xA; points_t[:,0,:] = in_particles.Points[contacts[:,1],:]&#xA; points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)&#xA; points.append(points_t)&#xA;&#xA;#merge everything&#xA;points = np.vstack(points)&#xA;vals = np.hstack(vals)&#xA;&#xA;#generate tubes&#xA;t = points[:,0,:]-points[:,1,:]&#xA;t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)&#xA;ez = np.where(t[:,1,None]&gt;t[:,2,None],np.array([[0,0,1]]),np.array([[0]]))&#xA;n1 = np.cross(t,ez)&#xA;n2 = np.cross(t,n1)&#xA;alphas = np.arange(0,2*np.pi, 2*np.pi/nf)&#xA;r = rf*vals**(1./2)*1&#xA;opoints = points[:,None,:,:] \&#xA; +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \&#xA; +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]&#xA;output.Points = opoints.reshape(-1,3)&#xA;output.PointData.append(np.repeat(vals,2*nf),&quot;Reaction&quot;)&#xA;n = points.shape[0]&#xA;pattern = np.ndarray([nf,4], np.int)&#xA;for i in range(nf) :&#xA; j = (i+1)%nf&#xA; pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)&#xA;types = np.full([n*nf],9,np.uint8)&#xA;locations = np.arange(0,5*n*nf,5,np.uint32)&#xA;cells = np.ndarray([n,nf,5],np.uint32)&#xA;cells[:,:,0] = 4&#xA;cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]&#xA;output.SetCells(types, locations, cells.reshape([-1]))&#xA;"
panel_visibility="advanced">
<Hints>
<Widget type="multi_line"/>
</Hints>
<Documentation>This property contains the text of a python program that
the programmable source runs.</Documentation>
</StringVectorProperty>
</SourceProxy>
</ProxyGroup>
</ServerManagerConfiguration>
......@@ -256,7 +256,7 @@ def getOutputDataSetTypeXml(info):
def getProxyGroup(info):
if not info.has_key("Group"):
if "Group" not in info:
return 'sources' if getNumberOfInputs(info) == 0 else 'filters'
else: return info["Group"]
......@@ -337,7 +337,7 @@ def replaceFunctionWithSourceString(namespace, functionName, allowEmpty=False):
def generatePythonFilterFromFiles(scriptFile, outputFile):
namespace = {}
execfile(scriptFile, namespace)
exec(compile(open(scriptFile).read(), scriptFile, 'exec'), namespace)
replaceFunctionWithSourceString(namespace, 'RequestData')
replaceFunctionWithSourceString(namespace, 'RequestInformation', allowEmpty=True)
......@@ -351,7 +351,7 @@ def generatePythonFilterFromFiles(scriptFile, outputFile):
def main():
if len(sys.argv) != 3:
print 'Usage: %s <python input filename> <xml output filename>' % sys.argv[0]
print(('Usage: %s <python input filename> <xml output filename>' % sys.argv[0]))
sys.exit(1)
inputScript = sys.argv[1]
......@@ -361,4 +361,4 @@ def main():
if __name__ == '__main__':
main()
\ No newline at end of file
main()
import numpy as np
nf = 5
rf = 2.5e-4
vals = inputs[0].CellData["Reaction"]
points = inputs[0].Points
cells = inputs[0].Cells.reshape([-1,3])[:,1:3]
mask = np.where(vals != 0)[0]
points = points[cells[mask,:],:]
vals = vals[mask]
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(t[:,1,None]>t[:,2,None],np.array([[0,0,1]]),np.array([[0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1
opoints = points[:,None,:,:] \
+n1[:,None,None,:]*r*np.sin(-alphas)[None,:,None,None] \
+n2[:,None,None,:]*r*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
j = (i+1)%nf
pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4