Commit 1ef2ef10 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

fix periodic

parent e06c3c58
Pipeline #8352 passed with stages
in 5 minutes and 3 seconds
......@@ -279,6 +279,16 @@ Mesh *mesh_load_msh(const char *filename)
m->parent_tag[nid] = m->parent_tag[nidSlave];
}
}
int *renumber = malloc(sizeof(int)*m->n_nodes);
int c = 0;
for (int i = 0; i < m->n_nodes; ++i) {
if (m->parent_tag[i] == i)
renumber[i] = c++;
}
for (int i = 0; i < m->n_nodes; ++i) {
m->parent_tag[i] = renumber[m->parent_tag[i]];
}
free(renumber);
free(boundary_tag_map);
mesh_gen_edges(m, n_bnd_el, bnd_elements, bnd_tags);
......
......@@ -341,25 +341,21 @@ class FluidProblem :
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()
rhs = np.zeros_like(solution_old)
self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
periodic_map = self._get_matrix("parent_tag", 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(self.elements(),self.n_fields(),self.solver_options
# )
sys = self.sys
rhs = np.zeros(sys.globalsize)
localv = np.ndarray([sys.localsize*self.n_elements()])
localm = np.ndarray([sys.localsize**2*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
sys.local_to_global(localv,localm,rhs)
if check_residual_norm > 0 :
norm0 = np.linalg.norm(rhs)
self.solution()[:] -= sys.solve(rhs)[periodic_map] # Matrice singuliere : logique car on mets plusieurs fois la meme ligne
# self.solution()[:] -= sys.solve(rhs)
self.solution()[:] -= sys.solve(rhs).reshape([-1,self.n_fields()])[periodic_map]
if check_residual_norm > 0 :
sys.zero_matrix()
rhs[:] = 0
......
......@@ -45,6 +45,7 @@ class LinearSystem :
self.localsize = elements.shape[1]*n_fields
self.elements = elements
self.idx = (self.elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
self.globalsize = nnodes*n_fields;
def local_to_global(self,localv,localm,rhs):
self.mat.zeroEntries()
......
......@@ -34,6 +34,8 @@ class LinearSystem :
self.cols = self.cols.reshape([-1])
self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
self.localsize = localsize
nnodes = np.max(elements)+1
self.globalsize = nnodes*n_fields;
def local_to_global(self,localv,localm,rhs):
np.add.at(rhs.reshape([-1]),self.idx,localv)
......
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