Commit 2d59104f authored by Michel Henry's avatar Michel Henry
Browse files

structure periodic_segments

parent d33f7f14
Pipeline #8404 passed with stages
in 3 minutes and 24 seconds
......@@ -921,3 +921,15 @@ size_t* gmsh_mesh_element_nodes_tag(GmshMesh *mesh, int dim, int ie, int ielt) {
}
return element;
}
int gmsh_mesh_n_periodic(GmshMesh *mesh) {
return mesh->n_periodic;
}
int gmsh_mesh_periodic_dim(GmshMesh *mesh, int ip) {
return mesh->periodic[ip]->dim;
}
int gmsh_mesh_periodic_tag(GmshMesh *mesh, int ip) {
return mesh->periodic[ip]->tag;
}
int gmsh_mesh_periodic_parent_tag(GmshMesh *mesh, int ip) {
return mesh->periodic[ip]->parent_tag;
}
......@@ -77,4 +77,9 @@ size_t gmsh_mesh_element_n_partitions(GmshMesh *mesh, int dim, int ie);
int* gmsh_mesh_element_partitions(GmshMesh *mesh, int dim, int ie);
size_t gmsh_mesh_element_n_nodes(GmshMesh *mesh, int dim, int ie);
size_t* gmsh_mesh_element_nodes_tag(GmshMesh *mesh, int dim, int ie, int ielt);
int gmsh_mesh_n_periodic(GmshMesh *mesh);
int gmsh_mesh_periodic_dim(GmshMesh *mesh, int ip);
int gmsh_mesh_periodic_tag(GmshMesh *mesh, int ip);
int gmsh_mesh_periodic_parent_tag(GmshMesh *mesh, int ip);
#endif
......@@ -423,7 +423,9 @@ class ParticleProblem :
def add_boundary_disk(self, x0, r, tag, material="default") :
self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
def add_boundary_periodic_segment(self, x0, x1, transform, tag) :
self._lib.particleProblemAddBoundaryPeriodicSegment.restype = c_size_t
return self._lib.particleProblemAddBoundaryPeriodicSegment(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*transform), tag.encode())
def add_boundary_segment(self, x0, x1, tag, material="default") :
self._lib.particleProblemAddBoundarySegment.restype = c_size_t
return self._lib.particleProblemAddBoundarySegment(self._p, self._coord_type(*x0), self._coord_type(*x1), tag.encode(),material.encode())
......@@ -514,6 +516,8 @@ class ParticleProblem :
phys_tag = self.physical_tags(_mesh, n_phys)
phys_dims = self.physical_dims(_mesh,n_phys)
physicals = [{},{},{},{}]
periodic_tag = self.periodic_tag(_mesh)
print(periodic_tag)
for i in range(n_phys):
physicals[int(phys_dims[i])][phys_names[i].decode()] = int(phys_tag[i])
entities = self.get_entities(_mesh)
......@@ -525,14 +529,23 @@ class ParticleProblem :
continue
tag = 0 if not ent.physicals else ent.physicals[0]
stag = self.getPhysicalName(physicals,ent.dimension, tag) or str(tag)
flag = tag in periodic_tag[ent.dimension]
print(f"tag = {tag}")
print(f"periodic_tag[{ent.dimension}] = {periodic_tag[ent.dimension]}")
print(f"\t\t => FLAG = {flag}")
for i, el in enumerate(ent.elements) :
for v in el:
v0 = el[1]
v1 = el[0]
if v[3] in self._addv :
continue
self._addv.add(v[3])
# if not flag : # influence du disque sur la frontiere periodique
self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
v0 = el[1]
v1 = el[0]
if flag:
transform = (0.0,0.0)
self.add_boundary_periodic_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), transform, stag)
continue
self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
else :
self._adde = set()
......@@ -563,6 +576,8 @@ class ParticleProblem :
(v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
(v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
stag,material)
if _mesh is not None :
self._lib.gmsh_mesh_free(_mesh)
def dim(self) :
return self._dim
......@@ -649,3 +664,19 @@ class ParticleProblem :
f = getattr(self._lib, "gmsh_mesh_element_nodes_tag")
f.restype = POINTER(c_size_t)
return None if nodes_by_elt < 1 else np.ctypeslib.as_array(f(_mesh, c_int(dim), c_int(ie), c_int(ielt)),shape=(nodes_by_elt,))
def n_periodic(self,_mesh):
self._lib.gmsh_mesh_n_periodic.restype = c_int
return self._lib.gmsh_mesh_n_periodic(_mesh)
def periodic_tag(self,_mesh):
n_periodic = self.n_periodic(_mesh)
vtag = [{},{},{},{}]
self._lib.gmsh_mesh_periodic_dim.restype = c_int
self._lib.gmsh_mesh_periodic_tag.restype = c_int
self._lib.gmsh_mesh_periodic_parent_tag.restype = c_int
for i in range(n_periodic):
dim = self._lib.gmsh_mesh_periodic_dim(_mesh,c_int(i))
tag = self._lib.gmsh_mesh_periodic_tag(_mesh,c_int(i))
parent_tag = self._lib.gmsh_mesh_periodic_parent_tag(_mesh,c_int(i))
vtag[int(dim)][int(tag)] = parent_tag
vtag[int(dim)][int(parent_tag)] = tag
return vtag
......@@ -2,26 +2,26 @@
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
*
*
* List of the contributors to the development of MigFlow: see AUTHORS file.
* Description and complete License: see LICENSE file.
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation, either version
* 3 of the License, or (at your option) any later version.
*
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program (see COPYING and COPYING.LESSER files). If not,
* along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>.
*/
#include "scontact.h"
#include "scontact.h"
#include "quadtree.h"
#include "vector.h"
#include <math.h>
......@@ -36,6 +36,7 @@ typedef struct _Particle Particle;
typedef struct _Contact Contact;
typedef struct _Disk Disk;
typedef struct _Segment Segment;
typedef struct _PeriodicSegment PeriodicSegment;
typedef struct _Triangle Triangle;
struct _ParticleProblem{
......@@ -48,6 +49,7 @@ struct _ParticleProblem{
int *particleMaterial;
int *ForcedFlag;
Segment *segments;
PeriodicSegment *periodicSegments;
#if FRICTION_ENABLED
#if ROTATION_ENABLED
double *omega;
......@@ -112,7 +114,7 @@ static void contactBuildBasis(Contact *c) {
#if DIMENSION==3
c->t[0] = -c->n[2];
c->t[1] = c->n[0];
c->t[2] = -c->n[1];
c->t[2] = -c->n[1];
double ddot = dot(c->n,c->t);
double norm = 0;
for(int i = 0;i<3;i++){
......@@ -228,12 +230,71 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
c->type = PARTICLE_DISK;
return c->D < alert;
}
struct _PeriodicSegment{
int tag;
double transform[DIMENSION]; //Pas le plus optimise
double p[2][DIMENSION];
};
static void periodicSegmentBoundingBox(const PeriodicSegment *s, double *pmin, double *pmax) {
for (int i = 0; i < DIMENSION; ++i) {
pmin[i] = fmin(s->p[0][i], s->p[1][i]);
pmax[i] = fmax(s->p[0][i], s->p[1][i]);
}
}
static int periodicSegmentProjectionIsInside(const void *ps, const double *x) {
const PeriodicSegment *s = (const PeriodicSegment*)ps;
double alpha = 0, beta = 0;
for (int i = 0; i < DIMENSION; ++i) {
const double d = s->p[0][i] - s->p[1][i];
alpha += (x[i] - s->p[0][i]) * d;
beta += (x[i] - s->p[1][i]) * d;
}
return (alpha < 0 && beta > 0);
}
static int periodicSegmentInitContact(size_t id, const PeriodicSegment *s, size_t particle, const Particle *p, double *x, double alert, Contact *c) {
c->o0 = id;
c->o1 = particle;
c->dv = 0;
c->ct = 0;
#if DIMENSION == 3
c->cs = 0;
#endif
double t[DIMENSION];
double nt2 = 0;
double dt = 0;
for (int i = 0; i <DIMENSION; ++i){
t[i] = s->p[1][i] - s->p[0][i];
dt += t[i] * (s->p[0][i] - x[i]);
nt2 += t[i] * t[i];
}
double nn2 = 0;
for (int i = 0; i < DIMENSION; ++i) {
c->n[i] = s->p[0][i] - x[i] - t[i] / nt2 * dt;
nn2 += c->n[i] * c->n[i];
}
const double nnorm = sqrt(nn2);
for (int i = 0; i < DIMENSION; ++i) {
c->n[i] /= nnorm;
}
contactBuildBasis(c);
c->D = nnorm - p->r;
if (c->D < 0 && periodicSegmentProjectionIsInside(s, x)) c->D = 0;
c->a0 = 0;
c->a1 = 1;
c->type = PARTICLE_SEGMENT;
return c->D >= 0 && c->D < alert;
}
struct _Segment{
Boundary b;
double p[2][DIMENSION];
};
static void coordAdd(double *x, double a, const double *y) {
x[0] += a * y[0];
x[1] += a * y[1];
......@@ -337,9 +398,9 @@ static int triangleProjectionIsInside(const void *pt, const double *x) {
double xp[3];
double dp[3] = {t->p[0][0] - x[0], t->p[0][1] - x[1], t->p[0][2] - x[2]};
const double nd = dot(n, dp);
for (int i = 0; i < 3; ++i)
for (int i = 0; i < 3; ++i)
xp[i] = x[i] + n[i] * nd;
double dx[3][3] =
double dx[3][3] =
{{xp[0] - t->p[0][0], xp[1] - t->p[0][1], xp[2] - t->p[0][2]},
{xp[0] - t->p[1][0], xp[1] - t->p[1][1], xp[2] - t->p[1][2]},
{xp[0] - t->p[2][0], xp[1] - t->p[2][1], xp[2] - t->p[2][2]}};
......@@ -402,7 +463,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
double ct = 0;
double vt = c->ct;
for(int i = 0;i<2;i++){
vt += (p->velocity[c->o0*2+i]-p->velocity[c->o1*2+i])*c->t[i];
vt += (p->velocity[c->o0*2+i]-p->velocity[c->o1*2+i])*c->t[i];
}
#if ROTATION_ENABLED
vt += p->omega[c->o0]*p->particles[c->o0].r+p->omega[c->o1]*p->particles[c->o1].r;
......@@ -460,7 +521,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1*2/7, c->t);
coordAdd(&p->velocity[c->o0 * DIMENSION], -cs*c->a0*2/7, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*c->a1*2/7, c->s);
for (int i = 0; i<3; i++){
for (int i = 0; i<3; i++){
p->omega[c->o0*3 + i] += (1/p->particles[c->o0].r)*c->a0*(cs*c->t[i]-ct*c->s[i])*5/7;
p->omega[c->o1*3 + i] += (1/p->particles[c->o1].r)*c->a1*(cs*c->t[i]-ct*c->s[i])*5/7;
}
......@@ -558,7 +619,7 @@ static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Bound
}
c->cs += cs;
#endif
#else
#else
ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct, c->t);
#if DIMENSION==3
......@@ -690,7 +751,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->ct = cold->ct;
#if DIMENSION == 3
c->cs = cold->cs;
#endif
#endif
#if ROTATION_ENABLED
#if DIMENSION==2
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0/3, c->t);
......@@ -808,7 +869,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
}
}
}
// Triangles
// Triangles
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
triangleBoundingBox(&p->triangles[i], amin, amax);
......@@ -824,7 +885,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->dv = cold->dv;
#if FRICTION_ENABLED
c->cs = cold->cs;
c->ct = cold->ct;
c->ct = cold->ct;
#if ROTATION_ENABLED
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*2/7, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*2/7, c->t);
......@@ -902,7 +963,7 @@ static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double t
#endif
while (! converged) {
#if FRICTION_ENABLED
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
return 0;
}
......@@ -1026,7 +1087,7 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
}
}
vectorClear(found);
#endif
#endif
if (close_to_bnd == 1) continue;
cellSearch(tree,amin,amax,&found);
for (size_t ifound = 0; ifound < vectorSize(found); ++ifound) {
......@@ -1056,7 +1117,7 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
#if DIMENSION == 3
fi += c->s[iD]*alpha*c->cs;
#endif
#endif
#endif
for (int jD = 0; jD < DIMENSION; jD++){
fj[iD*DIMENSION+jD] += fi*(-d*c->n[jD])/V;
}
......@@ -1086,7 +1147,7 @@ static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, dou
for (size_t ic = fifoPop(queue); ic != FIFO_EMPTY; ic=fifoPop(queue)){
Contact *c = &p->contacts[ic];
#if FRICTION_ENABLED
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
fifoFree(queue);
free(activeContact);
......@@ -1139,9 +1200,9 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
if(!p->ForcedFlag[j]){
oldVelocity[j*DIMENSION + i] = p->velocity[j * DIMENSION + i];
oldVelocity[j*DIMENSION + i] = p->velocity[j * DIMENSION + i];
p->velocity[j * DIMENSION + i] += externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
}
p->contactForces[j*DIMENSION+i] = p->velocity[j*DIMENSION+i];
}
......@@ -1208,7 +1269,15 @@ size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag, const char *material) {
return particleProblemAddBoundarySegmentTagId(p, x0, x1, particleProblemGetTagId(p, tag),particleProblemGetMaterialTagId(p, material));
}
size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double transform[DIMENSION], const char *tag){
PeriodicSegment *ps = vectorPush(&p->periodicSegments);
ps->tag = particleProblemGetTagId(p, tag);
for(int i = 0; i < DIMENSION; ++i){
ps->p[0][i] = x0[i];
ps->p[1][i] = x1[i];
ps->transform[i] = transform[i];
}
}
void particleProblemAddParticleMaterialTagId(ParticleProblem *p, const double x[DIMENSION], double r, double m, int tag, int forced) {
size_t n = vectorSize(p->particles);
Particle *particle = vectorPush(&p->particles);
......@@ -1256,7 +1325,7 @@ void particleProblemRemoveParticles(ParticleProblem *p, const int *keep_flag) {
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
if (keep_flag[i])
map[i] = nid++;
else
else
map[i] = -1;
}
int *keep_contacts = malloc(sizeof(int)*vectorSize(p->contacts));
......@@ -1386,6 +1455,7 @@ void particleProblemDelete(ParticleProblem *p) {
vectorFree(p->particles);
vectorFree(p->disks);
vectorFree(p->segments);
vectorFree(p->periodicSegments);
vectorFree(p->contacts);
vectorFree(p->savedContacts);
vectorFree(p->position);
......@@ -1466,6 +1536,7 @@ ParticleProblem *particleProblemNew() {
p->materialName = NULL;
p->disks = NULL;
p->segments = NULL;
p->periodicSegments = NULL;
p->particleMaterial = NULL;
p->position = NULL;
p->ForcedFlag = NULL;
......
......@@ -2,22 +2,22 @@
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
*
*
* List of the contributors to the development of MigFlow: see AUTHORS file.
* Description and complete License: see LICENSE file.
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation, either version
* 3 of the License, or (at your option) any later version.
*
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program (see COPYING and COPYING.LESSER files). If not,
* along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>.
*/
......@@ -36,6 +36,7 @@ int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, const char *material, int forced);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag, const char *material);
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag, const char *material);
size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double transform[DIMENSION], const char *tag);
size_t particleProblemGetTagId(ParticleProblem *p, const char *tag);
size_t particleProblemGetMaterialTagId(ParticleProblem *p, const char *tag);
const char* particleProblemGetTagName(ParticleProblem *p, size_t tag);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment