Commit 2b29a050 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

set_mesh_and_transfer_solution

parent bdda6b4f
Pipeline #5249 passed with stage
in 25 seconds
......@@ -1677,44 +1677,49 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
problem->boundaries = mesh_boundaries_create(problem->mesh);
}
void fluid_problem_adapt_load_mesh(FluidProblem *problem, const char *filename)
void fluid_problem_set_mesh_and_transfer_solution(FluidProblem *problem, Mesh *new_mesh)
{
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
new_eid[i] = -1;
double *newx = (double*)malloc(sizeof(double)*D*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < D; ++k)
newx[i*D+k] = new_mesh->x[i*3+k];
mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
for (int i = 0; i < new_mesh->n_nodes; ++i) {
unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
if(new_eid[i] < 0) {
printf("new mesh point not found in old mesh\n");
exit(1);
}
double phi[N_SF];
shape_functions(new_xi+i*D, phi);
for (int j=0; j<n_fields; ++j) {
new_solution[i*n_fields+j] = 0;
for (int k = 0; k < N_SF; ++k)
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
}
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
new_eid[i] = -1;
double *newx = (double*)malloc(sizeof(double)*D*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < D; ++k)
newx[i*D+k] = new_mesh->x[i*3+k];
mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
for (int i = 0; i < new_mesh->n_nodes; ++i) {
unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
if(new_eid[i] < 0) {
printf("new mesh point not found in old mesh\n");
exit(1);
}
free(new_eid);
free(new_xi);
free(newx);
fluid_problem_set_mesh(problem,new_mesh);
free(problem->solution);
problem->solution = new_solution;
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_porosity(problem);
double phi[N_SF];
shape_functions(new_xi+i*D, phi);
for (int j=0; j<n_fields; ++j) {
new_solution[i*n_fields+j] = 0;
for (int k = 0; k < N_SF; ++k)
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
}
}
free(new_eid);
free(new_xi);
free(newx);
fluid_problem_set_mesh(problem,new_mesh);
free(problem->solution);
problem->solution = new_solution;
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_porosity(problem);
}
void fluid_problem_adapt_load_mesh(FluidProblem *problem, const char *filename)
{
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
fluid_problem_set_mesh_and_transfer_solution(problem, new_mesh);
}
/*
......@@ -1755,7 +1760,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
problem->mesh = mesh;
}else if(boolean==1 && boolean_init==0){
printf("==> Call to the fluid_problem_adapt_gen_mesh function\n");
walk_tree(mesh);
//walk_tree(mesh);
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin, e_target, "adapt/lc.pos", bool_porosity, bool_pressure, bool_velocity);
fluid_problem_adapt_load_mesh(problem,"adapt/mesh.msh");
}
......@@ -1942,7 +1947,7 @@ double *fluid_problem_element_size(const FluidProblem *p)
return p->element_size;
}
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals)
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals, int transfer_solution)
{
Mesh *m = malloc(sizeof(Mesh));
m->n_elements = n_elements;
......@@ -1960,7 +1965,10 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
m->boundary_names[i] = strdup(physicals[i]);
}
mesh_gen_edges(m,n_boundaries, boundaries, boundary_tags);
fluid_problem_set_mesh(p,m);
if(transfer_solution)
fluid_problem_set_mesh_and_transfer_solution(p,m);
else
fluid_problem_set_mesh(p,m);
}
......
......@@ -124,7 +124,7 @@ int fluid_problem_n_nodes(const FluidProblem *p);
double *fluid_problem_coordinates(const FluidProblem *p);
int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals);
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals, int transfer_solution);
int fluid_problem_n_mesh_boundaries(const FluidProblem *p);
void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize);
void fluid_problem_mesh_boundary_interfaces(const FluidProblem *p, int bid, int *binterfaces);
......
......@@ -36,15 +36,6 @@ typedef struct {
char **boundary_names;
int n_interfaces;
int *interfaces; // eid0, closure0, {(eid1, closure1),(-1,phys_id)}
int n_phys;
char **phys_name;
int *phys_n_nodes;
int **phys_nodes;
int *phys_n_boundaries;
int **phys_boundaries; // eid, side
int *phys_dimension;
int *phys_id;
double *min_h_target; // added attribut for "in advance" strategy
} Mesh;
......
......@@ -288,7 +288,8 @@ class FluidProblem :
self._lib.fluid_problem_set_elements(self._fp,
c_int(x.shape[0]),np2c(x,np.float64),
c_int(el.shape[0]),np2c(el,np.int32),
c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names,
c_int(0)
)
sol = self.solution()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
......
......@@ -297,6 +297,7 @@ class ParticleProblem :
self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
_np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32))
mcheck = self._get_idx("ParticleMaterial")
self._fc = d["Contacts"][:,:self._dim]
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
......
......@@ -1211,8 +1211,8 @@ void particleProblemClearParticles(ParticleProblem *p) {
vectorClear(p->velocity);
#if FRICTION_ENABLED
vectorClear(p->omega);
vectorClear(p->particleMaterial);
#endif
vectorClear(p->particleMaterial);
}
void particleProblemAddParticles(ParticleProblem *p, int n_particles, double *position, double *mass, double *radius, double *velocity, double *omega, int *tag) {
......@@ -1241,9 +1241,15 @@ void particleProblemAddParticles(ParticleProblem *p, int n_particles, double *po
void particleProblemClearBoundaries(ParticleProblem *p) {
vectorClear(p->disks);
vectorClear(p->diskTag);
vectorClear(p->diskMaterial);
vectorClear(p->segments);
vectorClear(p->segmentTag);
vectorClear(p->segmentMaterial);
#if DIMENSION == 3
vectorClear(p->triangles);
vectorClear(p->triangleTag);
vectorClear(p->triangleMaterial);
#endif
}
......
......@@ -122,7 +122,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
......@@ -134,7 +134,7 @@ if init:
# Adaptation of the mesh.
#fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.adapt_gen_mesh(8*r, 4*r, 6602*2, "adapt/lc.pos")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
while t < tEnd1 :
# Fluid solver
......@@ -142,7 +142,7 @@ if init:
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# Increase time step if convergence
if niter < 5 and dt1 < 20:
dt1 *= 3
......@@ -154,7 +154,7 @@ if init:
t = 0
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#
# COMPUTATION LOOP
......@@ -188,7 +188,7 @@ while t < tEnd :
# load the save past state
fluid.import_vtk("output/fluid_00000.vtu")
p.read_vtk(outputdir,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# create the .msh on the adapt/lc.pos
fluid.adapt_load_mesh("adapt/mesh.msh")
#fluid.implicit_euler(dt)
......@@ -208,7 +208,7 @@ while t < tEnd :
p.iterate(dt/nsub, forces)
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# Output files writting
if ii %outf == 0 :
......
Markdown is supported
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