Commit 418d866e authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 7bec674d
Pipeline #9697 failed with stages
in 60 minutes and 40 seconds
......@@ -1187,8 +1187,15 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
const char *type_p[] = {"tetrahedron_p1",NULL};
#else
const char *types[] = {"triangle_p1","triangle_p1", "triangle_p1",NULL};
// const char *types[] = {"triangle_p2","triangle_p2", "triangle_p1",NULL};
const char *type_p[] = {"triangle_p1",NULL};
#endif
// const char *types[D+2];
// for(int i = 0; i < D+1; ++i){
// types[i] = *(((char**)etypes)+i);
// printf("type : %s\n", types[i]);
// }
// types[D+1] = NULL;
problem->fields = fe_fields_new(types);
problem->field_porosity = fe_fields_new(type_p);
const char *type_a[] = {"triangle_p1dg",NULL};
......@@ -1968,3 +1975,18 @@ void fluid_problem_reset_boundary_force(FluidProblem *problem) {
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
problem->boundary_force[i] = 0;
}
int fluid_problem_local_size(FluidProblem *problem) {
return problem->fields->local_size;
}
// int fluid_problem_global_size(FluidProblem *problem) {
// int n_nodes = problem->mesh->n_nodes;
// int n_edges = problem->mesh->n_edges;
// int n_faces = problem->mesh->n_faces;
// int n_elt = problem->mesh->n_elements;
// int global_size = problem->fields->nd[0]*n_nodes
// + problem->fields->nd[1]*n_edges
// + problem->fields->nd[2]*n_edges;
// #if DIMENSION == 3
// global_size += problem->fields->nd[3]*n_elt;
// #endif
// }
......@@ -111,7 +111,7 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor,int temporal, int advection);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection);
void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume);
void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh);
int *fluid_problem_particle_element_id(FluidProblem *problem);
......@@ -127,6 +127,7 @@ double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
int fluid_problem_local_size(FluidProblem *problem);
......
......@@ -89,7 +89,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_triangle_p1;
fe->df = fe_df_triangle_p1;
strcpy(fe->type, "triangle_p1");
fe->type = "triangle_p1";
}
else if (strcmp(etype,"triangle_p1dg") == 0) {
fe->nlocal = 3;
......@@ -99,7 +99,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_triangle_p1;
fe->df = fe_df_triangle_p1;
strcpy(fe->type, "triangle_p1dg");
fe->type = "triangle_p1dg";
}
else if (strcmp(etype,"triangle_p2") == 0) {
fe->nlocal = 6;
......@@ -109,7 +109,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_triangle_p2;
fe->df = fe_df_triangle_p2;
strcpy(fe->type, "triangle_p2");
fe->type = "triangle_p2";
}
else if (strcmp(etype,"line_p1") == 0) {
fe->nlocal = 2;
......@@ -119,7 +119,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_line_p1;
fe->df = fe_df_line_p1;
strcpy(fe->type, "line_p1");
fe->type = "line_p1";
}
else if (strcmp(etype,"line_p2") == 0) {
fe->nlocal = 3;
......@@ -129,7 +129,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_line_p2;
fe->df = fe_df_line_p2;
strcpy(fe->type, "line_p2");
fe->type = "line_p2";
}
else if (strcmp(etype,"tetrahedron_p1") == 0) {
fe->nlocal = 4;
......@@ -139,7 +139,7 @@ FEElement *fe_element_new(const char *etype) {
fe->n[3] = 0;
fe->f = fe_f_tetrahedron_p1;
fe->df = fe_df_tetrahedron_p1;
strcpy(fe->type, "tetrahedron_p1");
fe->type = "tetrahedron_p1";
}
else {
......
......@@ -4,13 +4,13 @@
#include <math.h>
#include "mesh.h"
#define deps 1
#define deps 1
struct FEElementStruct{
int nlocal;
int n[DIMENSION+1];
void (*f)(const double *xi, double *f);
void (*df)(const double *xi, const double dxidx[D][D], double f[][D]);
char type[256];
const char *type;
};
FEElement *fe_element_new(const char *etype);
......
......@@ -90,7 +90,7 @@ static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *
sort_edge_nodes(he->n);
he->element = -1;
he->pos = bnd_tags[i];
he ++;
he++;
}
if (periodic) {
for (int i = 0; i < n_half_edges; ++i){
......@@ -200,6 +200,8 @@ Mesh *mesh_new_from_elements(int n_nodes, double *x, int n_elements, int *elemen
for(int i = 0; i < m->n_nodes; ++i){
m->periodic_mapping[i] = periodic ? periodic[i] : i;
}
// add all edges
// add all faces
mesh_gen_edges(m, n_boundaries, boundaries, boundary_tags, 0, &m->n_interfaces, &m->interfaces);
mesh_gen_edges(m, n_boundaries, boundaries, boundary_tags, 1, &m->n_interfaces_periodic, &m->interfaces_periodic);
return m;
......
......@@ -133,7 +133,7 @@ def _load_msh(mesh_file_name, lib, dim) :
class FluidProblem :
"""Creates the numerical representation of the fluid."""
def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True):
def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True, etypes=None):
"""Builds the fluid structure.
Keyword arguments:
......@@ -164,6 +164,11 @@ class FluidProblem :
self._lib = lib3
else :
raise ValueError("dim should be 2 or 3.")
# if etypes is None and dim == 2:
# etypes = ["triangle_p1","triangle_p1", "triangle_p1"]
# if etypes is None and dim == 3:
# etypes = ["tetrahedron_p1","tetrahedron_p1","tetrahedron_p1","tetrahedron_p1"]
self._lib.fluid_problem_new.restype = c_void_p
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
......@@ -462,10 +467,10 @@ class FluidProblem :
solution_old = self.solution().copy()
self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
localsize = self._lib.fluid_problem_local_size(self._fp)
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
if self.sys is None :
self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options)
self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options, localsize)
sys = self.sys
rhs = np.zeros(sys.globalsize)
localv = np.ndarray([sys.localsize*self.n_elements()])
......
......@@ -24,7 +24,7 @@ from petsc4py import PETSc
import numpy as np
class LinearSystem :
def __init__(self, elements, n_fields, options) :
def __init__(self, elements, n_fields, options, localsize) :
nnodes = np.max(elements)+1
nn = elements.shape[1]
pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
......@@ -42,7 +42,7 @@ class LinearSystem :
self.ksp.setFromOptions()
self.mat.zeroEntries()
###
self.localsize = elements.shape[1]*n_fields
self.localsize = localsize
self.elements = elements
self.idx = (self.elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
self.globalsize = nnodes*n_fields;
......
......@@ -24,8 +24,7 @@ import scipy.sparse.linalg
class LinearSystem :
_initialized = False
def __init__(self, elements, n_fields,options) :
localsize = n_fields*elements.shape[1]
def __init__(self, elements, n_fields,options,localsize) :
self.el = elements
idx = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,localsize])
self.rows = np.repeat(idx[:,:,None],localsize,axis=2)
......
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