Skip to content

GitLab

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

fix mesh transfert

parent 253d240f
Pipeline #5314 passed with stage
in 25 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ a96dc4c0
......@@ -1796,10 +1796,13 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem){
int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int));
double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double));
int *boundaries_new = (int*) calloc(2*n_boundaries_new, sizeof(int));
int *boundary_tags_new = (int*) malloc(n_boundaries_new*sizeof(int));
// We keep a part of the vector x because no-nodes were removed (for the moment)
for (int ii=0; ii<mesh->n_nodes; ii++) {
x_new[ii*3+0] = mesh->x[ii*3+0]; x_new[ii*3+1] = mesh->x[ii*3+1];
x_new[ii*3+0] = mesh->x[ii*3+0];
x_new[ii*3+1] = mesh->x[ii*3+1];
x_new[ii*3+2] = 0.;
}
// initialize to zero
......@@ -1823,25 +1826,25 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem){
if (bool_boundary) {
boundaries_new[ibd*2+0] = node1;
boundaries_new[ibd*2+1] = node2;
boundary_tags_new[ibd] = mesh->interfaces[ii*4+3];
ibd++;
}
}
// Get the elements_new, x_new and boundaries_new arrays with the desired structure
create_the_new_mesh_with_interface(mesh, Tree, elements_new, n_elements_new, x_new, n_nodes_new, tag_node, boundaries_new, n_boundaries_new, ibd);
create_the_new_mesh_with_interface(mesh, Tree, elements_new, n_elements_new, x_new, n_nodes_new, tag_node, boundaries_new, boundary_tags_new,n_boundaries_new, ibd);
printf("HERE, creation of elements_new, x_new and boundaries_new\n");
// Creation of the new mesh and transfert the solution on the new mesh
int *boundary_tags = (int*) malloc(mesh->n_boundaries*sizeof(int));
for (int jj=0; jj<mesh->n_boundaries; jj++) {
boundary_tags[jj]=jj;
}
fluid_problem_set_elements(problem, n_nodes_new, x_new, n_elements_new, elements_new, n_boundaries_new, boundaries_new, boundary_tags, mesh->n_boundaries, mesh->boundary_names,1);
fluid_problem_set_elements(problem, n_nodes_new, x_new, n_elements_new, elements_new, n_boundaries_new, boundaries_new, boundary_tags_new, mesh->n_boundaries, mesh->boundary_names,1);
FILE *ftmp = fopen("test.msh","w");
mesh_write_msh(problem->mesh,ftmp);
fclose(ftmp);
free(Tree);
free(elements_new);
free(x_new);
free(boundaries_new);
free(boundary_tags);
free(boundary_tags_new);
}
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {
......@@ -2027,7 +2030,7 @@ double *fluid_problem_element_size(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, int transfer_solution)
{
printf("HERE, set mesh and transfert solution \n");
printf("HERE, set mesh and transfert solution %i elements\n", n_elements);
Mesh *m = malloc(sizeof(Mesh));
m->n_elements = n_elements;
m->elements = malloc(n_elements*sizeof(int)*N_N);
......@@ -2036,6 +2039,7 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
m->n_elements = n_elements;
m->x = malloc(n_nodes*3*sizeof(double));
memcpy(m->x,x,sizeof(double)*3*n_nodes);
printf("nodes :\n");
m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_elements);
memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);
m->n_boundaries = n_physicals;
......@@ -2043,6 +2047,10 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
for (int i = 0; i < n_physicals; ++i) {
m->boundary_names[i] = strdup(physicals[i]);
}
printf("n_boundaries : %i\n",n_boundaries);
for (int i = 0; i < n_boundaries; ++i) {
printf("%i %i with tag %i \n",boundaries[i*2+0],boundaries[i*2+1],boundary_tags[i]);
}
mesh_gen_edges(m,n_boundaries, boundaries, boundary_tags);
if(transfer_solution)
fluid_problem_set_mesh_and_transfer_solution(p,m);
......
fluid/mesh.c
View file @ a96dc4c0
......@@ -372,7 +372,18 @@ Mesh *mesh_load_msh(const char *filename)
int mesh_write_msh(Mesh *mesh, FILE *f)
{
int nbnd = 0;
for (int i = 0; i< mesh->n_interfaces; ++i) {
if (mesh->interfaces[i*4+2] == -1) nbnd++;
}
fprintf(f,"$MeshFormat\n2.2 0 0\n$EndMeshFormat\n");
fprintf(f,"$PhysicalNames\n");
fprintf(f,"%i\n", mesh->n_boundaries);
for (int i = 0; i < mesh->n_boundaries; ++i) {
fprintf(f,"%i %i \"%s\"\n", DIMENSION-1,i+1,mesh->boundary_names[i]);
}
fprintf(f,"$EndPhysicalNames\n");
fprintf(f,"$Nodes\n");
fprintf(f, "%d\n", mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i) {
......@@ -380,7 +391,7 @@ int mesh_write_msh(Mesh *mesh, FILE *f)
}
fprintf(f, "$EndNodes\n");
fprintf(f, "$Elements\n");
fprintf(f, "%i\n", mesh->n_elements);
fprintf(f, "%i\n", mesh->n_elements+nbnd);
for (int i= 0; i < mesh->n_elements; ++i) {
#if DIMENSION==2
int *tri = mesh->elements+i*3;
......@@ -390,6 +401,22 @@ int mesh_write_msh(Mesh *mesh, FILE *f)
fprintf(f, "%i 4 2 1 1 %i %i %i %i\n", i+1, tet[0]+1, tet[1]+1, tet[2]+1, tet[3]+1);
#endif
}
nbnd = 0;
for (int i = 0; i< mesh->n_interfaces; ++i) {
if (mesh->interfaces[i*4+2] != -1) continue;
int *cl = elbnd[mesh->interfaces[i*4+1]];
int iel = mesh->interfaces[i*4+0];
int phys = mesh->interfaces[i*4+3];
int imsh = nbnd+mesh->n_elements;
#if DIMENSION==2
int *tri = mesh->elements+iel*3;
fprintf(f, "%i 1 2 %i %i %i %i\n", imsh+1, phys+1,phys+1,tri[cl[0]]+1, tri[cl[1]]+1);
#else
int *tet = mesh->elements+iel*4;
fprintf(f, "%i 2 2 %i %i %i %i %i\n", imsh+1, phys+1,phys+1,tet[cl[0]]+1, tet[cl[1]]+1, tet[cl[2]]+1);
#endif
nbnd++;
}
fprintf(f, "$EndElements\n");
return 0;
}
......@@ -660,6 +687,7 @@ int need_refinement(double X1, double Y1, double X2, double Y2, double X3, doubl
/* This function returns the number of neighbours of the current element as well as their tag.
Remember that the tag of an element is its index in the array plus 1.
*/
/*
void walk_tree_with_interface(Mesh *mesh) {
// passes through the coarse mesh and create the Tree_cell structure
......@@ -684,7 +712,9 @@ void walk_tree_with_interface(Mesh *mesh) {
// We keep a part of the vector x because no-nodes were removed (for the moment)
for (int ii=0; ii<mesh->n_nodes; ii++) {
x_new[ii*3+0] = mesh->x[ii*3+0]; x_new[ii*3+1] = mesh->x[ii*3+1];
x_new[ii*3+0] = mesh->x[ii*3+0];
x_new[ii*3+1] = mesh->x[ii*3+1];
x_new[ii*3+2] = 0;
}
// initialize to zero
......@@ -713,12 +743,13 @@ void walk_tree_with_interface(Mesh *mesh) {
}
create_the_new_mesh_with_interface(mesh, Tree, elements_new, n_elements_new, x_new, n_nodes_new, tag_node, boundaries_new, n_boundaries_new, ibd);
free(Tree);
free(elements_new);
free(x_new);
free(boundaries_new);
}
*/
void create_the_tree_with_interface(Mesh *mesh, Tree_cell *Tree, int *n_elements_new, int *n_nodes_new, int *tag_node, int *n_boundaries_new){
......@@ -806,7 +837,7 @@ void create_the_tree_with_interface(Mesh *mesh, Tree_cell *Tree, int *n_elements
}
void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elements_new, int n_elements_new, double *x_new, int n_nodes_new, int tag_node, int *boundaries_new, int n_boundaries_new, int tag_boundary){
void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elements_new, int n_elements_new, double *x_new, int n_nodes_new, int tag_node, int *boundaries_new, int *boundary_tags_new, int n_boundaries_new, int tag_boundary){
// Initialized the index of nodes and elements
int node_index = tag_node;
......@@ -827,6 +858,7 @@ void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elemen
double X2 = mesh->x[node2*3]; double Y2 = mesh->x[node2*3+1];
x_new[Tree[iel].sub_nodes[ii]*3+0] = 0.5*(X1+X2);
x_new[Tree[iel].sub_nodes[ii]*3+1] = 0.5*(Y1+Y2);
x_new[Tree[iel].sub_nodes[ii]*3+2] = 0.;
}
}
}
......@@ -861,19 +893,21 @@ void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elemen
}else if (bool_one_node){
//printf("one node\n"); // information about boundaries will be found here
int tab[6] = {el[0], 0, el[1], 0, el[2], 0};
add_1_nodes(mesh, &Tree[iel], x_new, tab, XY, &node_index, boundaries_new, &tag_boundary, n_boundaries_new);
add_1_nodes(mesh, &Tree[iel], x_new, tab, XY, &node_index, boundaries_new, &tag_boundary, boundary_tags_new);
add_elements(tab, index, elements_new, &elt_index);
}else if(bool_two_nodes){
//printf("two nodes\n"); // information about boundaries will be found here
int tab[6] = {el[0], 0, el[1], 0, el[2], 0};
add_2_nodes(mesh, &Tree[iel], x_new, tab, XY, &node_index, boundaries_new, &tag_boundary);
add_2_nodes(mesh, &Tree[iel], x_new, tab, XY, &node_index, boundaries_new, boundary_tags_new, &tag_boundary);
add_elements(tab, index, elements_new, &elt_index);
}else{
//printf("three nodes\n");
for (int kk=0; kk<3; kk++) {
x_new[node_index*3+0] = XY[2*kk+1][0]; x_new[node_index*3+1] = XY[2*kk+1][1];
x_new[node_index*3+0] = XY[2*kk+1][0];
x_new[node_index*3+1] = XY[2*kk+1][1];
x_new[node_index*3+2] = 0;
node_index++;
}
int tab[6] = {el[0], node_index-4, el[1], node_index-3, el[2], node_index-2};
......@@ -891,6 +925,7 @@ void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elemen
for (int kk=0; kk<3; kk++) {
elements_new[elt_index*3+kk] = el[kk];
}
elt_index++;
}else if ((Tree[iel].sub_nodes[0]==-1 && Tree[iel].sub_nodes[1]==-1) || (Tree[iel].sub_nodes[1]==-1 && Tree[iel].sub_nodes[2]==-1) || (Tree[iel].sub_nodes[2]==-1 && Tree[iel].sub_nodes[0]==-1)){
//printf("Case2: the triangle should be cut into two pieces\n");
cut_into_2_with_interface(mesh, &Tree[iel], iel, elements_new, &elt_index);
......@@ -905,7 +940,24 @@ void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elemen
}
}
}
printing_post_processing_files(mesh, elements_new, x_new, boundaries_new, n_elements_new, n_boundaries_new);
//printing_post_processing_files(mesh, elements_new, x_new, boundaries_new, n_elements_new, n_boundaries_new);
printf("elt index : %i n_elements_new : %i\n", elt_index,n_elements_new);
// check orientation
for (int i= 0; i < n_elements_new; ++i) {
double *x[3];
for (int j = 0; j < 3; ++j) {
x[j] = &x_new[elements_new[i*3+j]*3];
}
double dx0[2] = {x[1][0]-x[0][0],x[1][1]-x[0][1]};
double dx1[2] = {x[2][0]-x[0][0],x[2][1]-x[0][1]};
double det = dx0[0]*dx1[1]-dx0[1]*dx1[0];
if (det < 0) {
int sw = elements_new[i*3];
elements_new[i*3] = elements_new[i*3+1];
elements_new[i*3+1] = sw;
}
}
}
void add_elements(int tab[6], int index[4][3], int *elements_new, int *elt_index){
......@@ -917,6 +969,21 @@ void add_elements(int tab[6], int index[4][3], int *elements_new, int *elt_index
}
}
void add_boundary_if_needed(int *boundaries_new, int *boundary_tags_new, int *tag_boundary, int node0, int node1, int midnode)
{
// Add the edges that appear to be on the boundary
int ibd = find_tag_in_boundaries(boundaries_new, node0, node1, *tag_boundary);
if (ibd == -1) return;
boundaries_new[ibd*2+0] = node0;
boundaries_new[ibd*2+1] = midnode;
boundaries_new[*tag_boundary*2+0] = midnode;
boundaries_new[*tag_boundary*2+1] = node1;
boundary_tags_new[*tag_boundary] = boundary_tags_new[ibd];
*tag_boundary+=1;
}
/*
This function add one node to the array x_new and create the array tab[] with the connection between the nodes and the created ones for elements construction.
This function is typically used if the elements has two boundary edges.
......@@ -926,7 +993,7 @@ void add_elements(int tab[6], int index[4][3], int *elements_new, int *elt_index
int tab[6], the node tag to create the elements
int *node_index, argument passed by reference
*/
void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary, int n_boundaries_new){
void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary, int *boundary_tags_new){
for (int ii=0; ii<3; ii++) {
if (current->sub_nodes[ii]!=-1 && current->sub_nodes[(ii+1)%3]!=-1) {
// Description of the first edge
......@@ -959,6 +1026,7 @@ void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], doub
int tag3 = *node_index;
x_new[*node_index*3+0] = 0.5*(X2e1+X2e2)*bool_1e1_1e2 + 0.5*(X2e1+X1e2)*bool_1e1_2e2 + 0.5*(X1e1+X2e2)*bool_2e1_1e2 + 0.5*(X1e1+X1e2)*bool_2e1_2e2;
x_new[*node_index*3+1] = 0.5*(Y2e1+Y2e2)*bool_1e1_1e2 + 0.5*(Y2e1+Y1e2)*bool_1e1_2e2 + 0.5*(Y1e1+Y2e2)*bool_2e1_1e2 + 0.5*(Y1e1+Y1e2)*bool_2e1_2e2;
x_new[*node_index*3+2] = 0.;
// Create the tab[] table
tab[0] = (int) (0.5*(node1e1+node1e2)*bool_1e1_1e2 + 0.5*(node1e1+node2e2)*bool_1e1_2e2 + 0.5*(node2e1+node1e2)*bool_2e1_1e2 + 0.5*(node2e1+node2e2)*bool_2e1_2e2);
......@@ -968,15 +1036,7 @@ void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], doub
tab[4] = node2e2*bool_1e1_1e2 + node1e2*bool_1e1_2e2 + node2e2*bool_2e1_1e2 + node1e2*bool_2e1_2e2;
tab[5] = tag2;
*node_index+=1;
// Add the edges that appear to be on the boundary
int ibd = find_tag_in_boundaries(boundaries_new, tab[2], tab[4], n_boundaries_new);
boundaries_new[ibd*2+0] = tab[2];
boundaries_new[ibd*2+1] = tag3;
boundaries_new[*tag_boundary*2+0] = tag3;
boundaries_new[*tag_boundary*2+1] = tab[4];
*tag_boundary+=1;
add_boundary_if_needed(boundaries_new, boundary_tags_new, tag_boundary, tab[2],tab[4],tag3);
}
}
}
......@@ -999,7 +1059,7 @@ int find_tag_in_boundaries(int *boundaries_new, int tag1, int tag2, int n_bounda
int tab[6], the node tag to create the elements
int *node_index, argument passed by reference
*/
void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary){
void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *boundary_tags_new, int *tag_boundary){
for (int ii=0; ii<3; ii++) {
if (current->sub_nodes[ii]==-1 && current->sub_nodes[(ii+1)%3]==-1) {
tab[1] = current->sub_nodes[(ii+2)%3];
......@@ -1007,7 +1067,7 @@ void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], doub
int tri_ind = mesh->interfaces[edge_ind*4+0];
int closure = mesh->interfaces[edge_ind*4+1];
const unsigned int *buf = &mesh->elements[tri_ind*N_N];
const int *buf = &mesh->elements[tri_ind*N_N];
int node1 = buf[elbnd[closure][0]];
int node2 = buf[elbnd[closure][1]];
......@@ -1018,16 +1078,28 @@ void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], doub
int bool_edge01 = ((bool01&&bool12) || (bool02&&bool11));
int bool_edge12 = ((bool11&&bool22) || (bool12&&bool12));
int bool_edge20 = ((bool21&&bool02) || (bool22&&bool01));
printf("%i %i %i\n", bool_edge01, bool_edge12, bool_edge20);
{
tab[3] = *node_index;
x_new[*node_index*3+0] = XY[1][0]*bool_edge20 + XY[3][0]*bool_edge01 + XY[5][0]*bool_edge12;
x_new[*node_index*3+1] = XY[1][1]*bool_edge20 + XY[3][1]*bool_edge01 + XY[5][1]*bool_edge12;
x_new[*node_index*3+2] = 0.;
*node_index+=1;
int node0 = tab[0]*bool_edge20 + tab[2]*bool_edge01 + tab[4]*bool_edge12;
int node1 = tab[2]*bool_edge20 + tab[4]*bool_edge01 + tab[0]*bool_edge12;
add_boundary_if_needed(boundaries_new, boundary_tags_new, tag_boundary, node0,node1,tab[3]);
}
tab[3] = *node_index;
x_new[*node_index*3+0] = XY[1][0]*bool_edge20 + XY[3][0]*bool_edge01 + XY[5][0]*bool_edge12;
x_new[*node_index*3+1] = XY[1][1]*bool_edge20 + XY[3][1]*bool_edge01 + XY[5][1]*bool_edge12;
*node_index+=1;
tab[5] = *node_index;
x_new[*node_index*3+0] = XY[1][0]*bool_edge12 + XY[3][0]*bool_edge20 + XY[5][0]*bool_edge01;
x_new[*node_index*3+1] = XY[1][1]*bool_edge12 + XY[3][1]*bool_edge20 + XY[5][1]*bool_edge01;
*node_index+=1;
{
tab[5] = *node_index;
x_new[*node_index*3+0] = XY[1][0]*bool_edge12 + XY[3][0]*bool_edge20 + XY[5][0]*bool_edge01;
x_new[*node_index*3+1] = XY[1][1]*bool_edge12 + XY[3][1]*bool_edge20 + XY[5][1]*bool_edge01;
x_new[*node_index*3+2] = 0;
*node_index+=1;
int node0 = tab[4]*bool_edge20 + tab[4]*bool_edge01 + tab[0]*bool_edge12;
int node1 = tab[2]*bool_edge20 + tab[0]*bool_edge01 + tab[2]*bool_edge12;
add_boundary_if_needed(boundaries_new, boundary_tags_new, tag_boundary, node0,node1,tab[5]);
}
}
}
}
......@@ -1221,6 +1293,11 @@ void cut_into_4_with_interface(Mesh *mesh, Tree_cell *current, int iel, int *ele
elements_new[*elt_index*3+1] = node1e3*(bool_1e3_1e1+bool_1e3_2e1)+node2e3*(bool_2e3_1e1+bool_2e3_2e1);
elements_new[*elt_index*3+2] = tag12;
*elt_index+=1;
elements_new[*elt_index*3+0] = tag12;
elements_new[*elt_index*3+1] = tag23;
elements_new[*elt_index*3+2] = tag31;
*elt_index+=1;
}
void printing_post_processing_files(Mesh *mesh, int *elements_new, double *x_new, int *boundaries_new, int n_elements_new, int n_boundaries_new){
......
fluid/mesh.h
View file @ a96dc4c0
......@@ -152,11 +152,11 @@ void cut_into_4_with_interface(Mesh *mesh, Tree_cell *current, int iel, int *ele
int find_tag_in_boundaries(int *boundaries_new, int tag1, int tag2, int n_boundaries_new);
void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elements_new, int n_elements_new, double *x_new, int n_nodes_new, int tag_node, int *boundaries_new, int n_boundaries_new, int tag_boundary);
void create_the_new_mesh_with_interface(Mesh *mesh, Tree_cell *Tree, int *elements_new, int n_elements_new, double *x_new, int n_nodes_new, int tag_node, int *boundaries_new, int *boundary_tags_new, int n_boundaries_new, int tag_boundary);
void create_the_tree_with_interface(Mesh *mesh, Tree_cell *Tree, int *n_elements_new, int *n_nodes_new, int *tag_node, int *n_boundaries_new);
void add_elements(int tab[6], int index[4][3], int *elements_new, int *elt_index);
void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary, int n_boundaries_new);
void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary);
void add_1_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *tag_boundary, int *boundary_tags_new);
void add_2_nodes(Mesh *mesh, Tree_cell *current, double *x_new, int tab[6], double XY[6][2], int *node_index, int *boundaries_new, int *boundary_tags_new, int *tag_boundary);
void printing_post_processing_files(Mesh *mesh, int *elements_new, double *x_new, int *boundaries_new, int n_elements_new, int n_boundaries_new);
#endif
......
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