Commit 56557cd2 authored by Michel Henry's avatar Michel Henry
Browse files

generalise mapping

parent 296c30bd
Pipeline #9958 failed with stages
in 2 minutes and 43 seconds
......@@ -478,7 +478,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
const Mesh *mesh = problem->mesh;
const size_t local_size = problem->fields->local_size;
const FEFields *fields = problem->fields;
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
......@@ -494,9 +494,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
int n_value = weak_boundary_n_values(wbnd);
double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
weak_boundary_values(problem->fields, mesh, bnd, wbnd, data);
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
const int iel = interface[0];
double *local_vector = &all_local_vector[local_size*iel];
......@@ -1138,6 +1136,7 @@ void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
}
}
// void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
// {
// const Mesh *mesh = problem->mesh;
......@@ -2012,7 +2011,7 @@ void weak_boundary_values(const FEFields *fields, const Mesh *mesh, const MeshBo
for(int id=0; id < D; ++id){
double xx = 0;
for(int isf=0; isf < mesh_element->nlocal; ++isf){
xx += mesh->x[el[isf]*3+id]*sf[isf];
xx += mesh->x[el[isf]*mesh_element->nlocal+id]*sf[isf];
}
x[(i*N_LQP+iqp)*D+id] = xx;
}
......@@ -2031,21 +2030,25 @@ 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_faces;
// #if DIMENSION == 3
// global_size += problem->fields->nd[3]*n_elt;
// #endif
// }
void fluid_problem_fields_local_map(FluidProblem *p, int *idx) {
p->fields->map = idx;
}
void fluid_problem_concentration_local_map(FluidProblem *p, int *idx) {
p->field_concentration->map = idx;
}
void fluid_problem_porosity_local_map(FluidProblem *p, int *idx) {
p->field_porosity->map = idx;
}
void fluid_problem_mesh_local_map(FluidProblem *p, int *idx) {
p->mesh->coord_fields->map = idx;
}
void fluid_problem_local_map(FluidProblem *p, int map[]){
for(int iel = 0; iel < p->mesh->n_elements; ++iel){
fe_fields_local_map(p->fields, p->mesh, iel, &(map[iel*p->fields->local_size]));
void fluid_problem_local_map(FluidProblem *p, int *idx) {
for (int iel = 0; iel < p->mesh->n_elements; ++iel){
int localsize = p->fields->local_size;
fe_fields_local_map(p->fields, p->mesh, iel, &(idx[iel*localsize]));
}
}
\ No newline at end of file
......@@ -127,8 +127,12 @@ 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);
void fluid_problem_local_map(FluidProblem *p, int map[]);
int fluid_problem_local_size(FluidProblem *p);
void fluid_problem_fields_local_map(FluidProblem *p, int *idx);
void fluid_problem_mesh_local_map(FluidProblem *p, int *idx);
void fluid_problem_concentration_local_map(FluidProblem *p, int *idx);
void fluid_problem_porosity_local_map(FluidProblem *p, int *idx);
void fluid_problem_local_map(FluidProblem *p, int *idx);
// mesh data (for i/o)
......
......@@ -182,6 +182,7 @@ FEFields *fe_fields_new(const char **etypes) {
FEFields *fields = malloc(sizeof(FEFields));
fields->n = 0;
fields->local_size = 0;
fields->map = NULL;
for (const char **etype = etypes; *etype != NULL; etype++) {
fields->n++;
}
......@@ -224,29 +225,41 @@ void fe_fields_df(const FEFields *fields, const double *xi, const double dxidx[D
}
}
void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]) {
void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]) { // unused
int c[fields->n];
c[0] = 0;
for (int i = 1; i < fields->n; ++i) {
c[i] = c[i-1]+fields->element[i-1]->nlocal;
}
int l = 0;
int n_n = DIMENSION+1;
for (int i = 0; i < fields->n; ++i) {
for (int j = 0; j < fields->element[i]->n[0]; ++j){
for (int k = 0; k < DIMENSION+1; ++k) {
map[c[i]++] = mesh->elements[iel*(DIMENSION+1)+k]*fields->nd[0]+l;
for (int k = 0; k < n_n; ++k) {
map[c[i]++] = mesh->elements[iel*n_n+k]*fields->nd[0]+l;
}
l++;
}
}
int shift = l*mesh->n_nodes;
// #if DIMENSION == 2
// int n_e = 3;
// #else
// int n_e = 6;
// #endif
// int shift = l*mesh->n_nodes;
// int edges[n_e *mesh->n_elements];
// mesh_gen_edges_by_element(mesh, edges);
// todo add edges (and faces in 3D)
for (int i = 0; i < fields->n; ++i) {
for (int j = 0; j < fields->element[i]->n[DIMENSION]; ++j){
map[c[i]++] = shift+iel*fields->nd[DIMENSION]+l;
l++;
}
}
// l = 0;
// for (int i = 0; i < fields->n; ++i){
// for(int j = 0; j < fields->element[i]->n[1]; ++j){
// for(int k = 0; k < n_e; ++k){
// map[c[i]++] = shift+edges[iel*n_e+k]*fields->nd[1]+l;
// }
// l++;
// }
// }
// unused
if (fields->n == 3) {
int *el = mesh->elements + 3*iel;
......@@ -259,8 +272,7 @@ void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int
}
void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]) {
int map[fields->local_size];
fe_fields_local_map(fields, mesh, iel, map);
int *map = &(fields->map[iel*fields->local_size]);
int c = 0;
for (int i = 0; i < fields->n; ++i) {
s[i] = 0;
......@@ -272,8 +284,7 @@ void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const
}
void fe_fields_eval_grad_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double dsf[][DIMENSION], const double *v, double s[], double ds[][DIMENSION]) {
int map[fields->local_size];
fe_fields_local_map(fields, mesh, iel, map);
int *map = &(fields->map[iel*fields->local_size]);
int c = 0;
for (int i = 0; i < fields->n; ++i) {
s[i] = 0;
......
......@@ -22,6 +22,7 @@ struct FEFieldsStruct{
int local_size;
int n;
int nd[DIMENSION+1];
int *map;
FEElement **element;
};
......
......@@ -199,15 +199,12 @@ 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;
}
MeshBoundary *mesh_boundaries_create(const Mesh *m) {
// To extend for p2p1
MeshBoundary *b = malloc(sizeof(MeshBoundary)*m->n_boundaries);
for (int i = 0; i < m->n_boundaries; ++i) {
b[i].n_nodes = 0;
......
......@@ -49,7 +49,6 @@ typedef struct {
void mesh_free(Mesh *m);
Mesh *mesh_new_from_elements(int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals, const int *periodic);
// not really necessary but useful to implement bnd conditions,
// can be rebuilt from the interface list
typedef struct {
......
......@@ -31,7 +31,7 @@ set(SRC
petsc4pylsys.py
scipylsys.py
VTK.py
# dof_numbering.py
dof_numbering.py
)
if(BUILD_MUMPS)
......
import numpy as np
import gmsh
gmsh.initialize()
from ._tools import gmsh
# gmsh.initialize()
#static p1 arrays, could be generated automatically but hard-coded to match gmsh order
tri_nodes = np.array([[0],[1],[2]])
......
......@@ -33,8 +33,11 @@ from time import time
import numpy as np
import os
import sys
from . import VTK
from ._tools import gmsh, get_linear_system_package, _np2c, timeit
from .dof_numbering import MeshCl, LagrangeSpace
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
......@@ -166,6 +169,7 @@ class FluidProblem :
if self._fp == None :
raise NameError("Cannot create fluid problem.")
def __del__(self):
"""Deletes the fluid structure."""
if(self._fp is not None) :
......@@ -182,6 +186,8 @@ class FluidProblem :
"""
mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
self._lib.fluid_problem_set_mesh(self._fp, mesh)
idx = self.set_numbering([1,1])
self._lib.fluid_problem_mesh_local_map(self._fp, _np2c(idx, dtype=np.int32))
self.sys = None
gmsh.model.remove()
......@@ -427,14 +433,20 @@ class FluidProblem :
self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
self._lib.fluid_problem_apply_boundary_conditions(self._fp)
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))
idx = self.set_numbering([1,1,1])
idx_p1 = self.set_numbering([1])
self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_p1, dtype=np.int32))
self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_p1, dtype=np.int32))
# print(time()-timer)
idx = np.zeros((self.n_elements()*self.local_size()), dtype=np.int32)
self.local_map(idx)
# idx = np.zeros((self.n_elements()*self.local_size()), dtype=np.int32)
# self.local_map(idx)
idx = idx.reshape([self.n_elements(), -1])
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
solution_old = self.solution().copy()
if self.sys is None :
constraints = []
if self.mean_pressure is not None:
......@@ -638,9 +650,35 @@ class FluidProblem :
self._lib.fluid_problem_local_size.restype = c_int
return self._lib.fluid_problem_local_size(self._fp)
def local_map(self, idx):
self._lib.fluid_problem_local_map(self._fp, _np2c(idx, dtype=np.int32))
@timeit
def local_map(self, map):
self._lib.fluid_problem_local_map(self._fp, _np2c(map, dtype=np.int32))
def set_numbering(self, order=[]):
mesh_cl = MeshCl(self.elements())
maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
ndofs = [lmap.shape[1] for lmap in maps]
ndofs_unique = np.unique(ndofs)
n_order = len(order)
# map en u0,v0,p0, u1,v1,p1, ...
idx_uvp = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
for i,ndof in enumerate(ndofs_unique):
ndof_by_order = ndof if i == 0 else ndof-ndofs_unique[i-1]
offset_dof = 0 if i == 0 else idx_uvp.max()+1
offset_ind = np.where(idx_uvp[0]==-1)[0][0]
order_to_use = np.where(ndofs >= ndof)[0]
offset_order = 0 if i == 0 else ndof-ndofs_unique[i-1]
n_order = len(order_to_use)
for j in range(ndof_by_order):
for k in range(n_order):
lmap = maps[order_to_use[k]][:,offset_order:] - maps[order_to_use[k]][:,offset_order:].min()
idx_uvp[:,offset_ind + j*n_order+k] = n_order*lmap[:,j]+k+offset_dof
# map en u0,u1, ..., v0, v1, ..., p0, p1, ...
idx = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
n = ndofs_unique[0]
for i in range(n_order):
idx[:,i*n:(i+1)*n] = idx_uvp[:,i::n_order]
return idx.reshape([-1])
def _n_physicals(self):
f = self._lib.fluid_problem_private_n_physical
......
......@@ -73,10 +73,6 @@ class Poiseuille(unittest.TestCase) :
fluid.set_open_boundary("RightUp",pressure=0)
fluid.set_open_boundary("RightDown",pressure=0)
fluid.set_strong_boundary("LeftUp", 1, 0)
fluid.set_strong_boundary("LeftUp", 0, lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftDown", 0, lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftDown", 1, 0)
ii = 0
t = 0
......@@ -91,7 +87,7 @@ class Poiseuille(unittest.TestCase) :
tic = time.perf_counter()
while ii < 100 :
#Fluid solver
time_integration.iterate(fluid,None,dt,check_residual_norm=-1)
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-5)
t += dt
#Output files writting
if ii %outf == 0 :
......
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