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

broken

parent ce0e915c
Pipeline #5505 passed with stage
in 2 minutes and 20 seconds
......@@ -22,6 +22,7 @@
*/
#include "tools.h"
#include "vector.h"
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
......@@ -1147,96 +1148,17 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
}
}
#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.02);
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->mesh->n_nodes; ++i){
problem->porosity[i] = 0.;
}
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;
int nclose = 0;
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] >= 0;
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;
nclose ++;
}
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 (nclose >= 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 j = problem->particle_nodes_n[i]; j < problem->particle_nodes_n[i+1]; ++j) {
problem->porosity[problem->particle_nodes[j]] += problem->particle_volume[i]*problem->particle_nodes_weight[j];
}
}
for (int i = 0; i < mesh->n_nodes; ++i){
for (int i = 0; i < problem->mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
problem->porosity[i] = 0.;
}
......@@ -1244,10 +1166,6 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
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)
......@@ -1303,6 +1221,9 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
problem->particle_contact = NULL;
problem->volume_drag = volume_drag;
problem->boundaries = NULL;
problem->particle_nodes_weight = NULL;
problem->particle_nodes = NULL;
problem->particle_nodes_n = NULL;
return problem;
}
......@@ -1316,6 +1237,7 @@ void fluid_problem_set_reduced_gravity(FluidProblem *p, int reduced_gravity) {
void fluid_problem_free(FluidProblem *problem)
{
free(problem->particle_nodes_n);
free(problem->mu);
free(problem->rho);
free(problem->solution);
......@@ -1344,6 +1266,8 @@ void fluid_problem_free(FluidProblem *problem)
mesh_tree_free(problem->mesh_tree);
mesh_free(problem->mesh);
mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
vectorFree(problem->particle_nodes);
vectorFree(problem->particle_nodes_weight);
free(problem);
}
......@@ -1708,6 +1632,103 @@ void fluid_problem_after_import(FluidProblem *problem)
#endif
}
static void compute_particle_nodes_weights(FluidProblem *problem) {
Mesh *mesh = problem->mesh;
int *touched_el_flag = malloc(sizeof(int)*mesh->n_elements);
double L2 = pow2(0.02);
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){
nodes_w[i] = -1;
}
vectorClear(problem->particle_nodes);
vectorClear(problem->particle_nodes_weight);
for (int i = 0; i < problem->n_particles; ++i) {
problem->particle_nodes_n[i] = vectorSize(problem->particle_nodes);
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;
int nclose = 0;
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] >= 0;
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;
nclose ++;
}
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 (nclose >= 3) {
for (size_t j = 0; j < vectorSize(nodes); ++j) {
*vectorPush(&problem->particle_nodes_weight) = nodes_w[nodes[j]]/wtot;
*vectorPush(&problem->particle_nodes) = nodes[j];
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){
*vectorPush(&problem->particle_nodes_weight) = sf[j];
*vectorPush(&problem->particle_nodes) = el[j];
}
}
for (int j = 0; j < vectorSize(stack); ++j) {
touched_el_flag[stack[j]] = 0;
}
}
problem->particle_nodes_n[problem->n_particles] = vectorSize(problem->particle_nodes);
free(neighbours);
vectorFree(nodes);
vectorFree(stack);
free(nodes_w);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload)
{
......@@ -1718,6 +1739,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
if(problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
free(problem->particle_nodes_n);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
......@@ -1725,18 +1747,18 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
free(problem->particle_volume);
free(problem->particle_velocity);
free(problem->particle_contact);
problem->particle_uvw = (double*)malloc(sizeof(double)*D*n);
problem->particle_element_id = (int*)malloc(sizeof(int)*n);
problem->particle_uvw = malloc(sizeof(double)*D*n);
problem->particle_element_id = malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=-1;
}
problem->particle_position = (double*)malloc(sizeof(double)*D*n);
problem->particle_mass = (double*)malloc(sizeof(double)*n);
problem->particle_force = (double*)malloc(sizeof(double)*n*D);
printf("alloc particle force : %i\n", n*D);
problem->particle_volume = (double*)malloc(sizeof(double)*n);
problem->particle_velocity = (double*)malloc(sizeof(double)*n*D);
problem->particle_contact = (double*)malloc(sizeof(double)*n*D);
problem->particle_position = malloc(sizeof(double)*D*n);
problem->particle_mass = malloc(sizeof(double)*n);
problem->particle_force = malloc(sizeof(double)*n*D);
problem->particle_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*D);
problem->particle_contact = malloc(sizeof(double)*n*D);
problem->particle_nodes_n = malloc(sizeof(int)*(n+1));
}
if (elid) {
for (int i = 0; i < n; ++i) {
......@@ -1758,6 +1780,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->oldporosity[i] = problem->porosity[i];
}
}
compute_particle_nodes_weights(problem);
fluid_problem_compute_porosity(problem);
if (init){
for (int i=0; i<problem->mesh->n_nodes; ++i){
......
......@@ -85,6 +85,11 @@ struct FluidProblem {
double *particle_contact;
double *particle_mass;
double *particle_force;
double *particle_node_forces;
double *particle_node_forces_du;
double *particle_nodes_weight;
int *particle_nodes;
int *particle_nodes_n;
int n_fluids;
//stabilisation coefficients
double *taup;
......
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