Commit 31612c99 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

weak_bnd : wip

parent 27b7b0db
Pipeline #4728 passed with stage
in 29 seconds
......@@ -316,7 +316,7 @@ static void f_wall(FluidProblem *problem,const double *n, double *f,const double
f[P] = 0;
for (int id = 0; id < D; ++id) {
f[U+id] = (v==NULL ? p : v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
f[U+id] = (v==NULL?p:v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
}
}
......@@ -443,53 +443,48 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
double *g0 = malloc(sizeof(double)*n_fields);
double *h0 = malloc(sizeof(double)*n_fields*D);
int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
f_cb *f;
#if DIMENSION == 2
int elbnd[3][2] = {{0,1},{1,2},{2,0}};
#else
int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
for (int i =0; i < mesh->n_nodes; ++i)
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
for (int i = 0; i < mesh->n_boundaries; ++i) {
if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
bndid = i;
}
}
if (bndid == -1) {
printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
}
MeshBoundary *bnd = &problem->boundaries[bndid];
double *v = malloc(wbnd->n_value*bnd->n_nodes*sizeof(double));
for (int i = 0; i < mesh->n_nodes; ++i) {
bnd_node_local_id[i] = -1;
if(mesh->phys_dimension[iphys] != D-1) continue;
double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
int j = mesh->phys_nodes[iphys][i];
bnd_node_local_id[j] = i;
for (int k = 0; k < D; ++k)
x[i*D+k] = mesh->x[j*3+k];
}
f = NULL;
int n_value;
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *bnd = problem->weak_boundaries + ibnd;
if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
f = bnd->cb;
n_value = bnd->n_value;
if (n_value > 0){
v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
}
break;
if (wbnd->n_value != 0) {
double *x = (double*)malloc(bnd->n_nodes*sizeof(double)*D);
for (int i = 0; i < bnd->n_nodes; ++i){
int j = bnd->nodes[i];
bnd_node_local_id[j] = i;
for (int k = 0; k < D; ++k)
x[i*D+k] = mesh->x[j*3+k];
}
wbnd->apply(bnd->n_nodes, x, v);
free(x);
}
if (f == NULL) {
printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
exit(1);
}
for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
const int iel = bnd[0];
const int iebnd = bnd[1];
double *vbnd = malloc(sizeof(double)*wbnd->n_value);
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
const int iel = interface[0];
const int icl = interface[1];
int nodes[D];
double *x[D];
const unsigned int *el = &mesh->elements[iel*N_N];
for (int i = 0; i < D; ++i){
nodes[i] = el[elbnd[iebnd][i]];
nodes[i] = el[elbnd[icl][i]];
x[i] = &mesh->x[nodes[i]*3];
}
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
......@@ -516,7 +511,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
double xi[D];
#if DIMENSION == 2
double a = LQP[iqp];
switch(iebnd) {
switch(icl) {
case 0 : xi[0] = a; xi[1] = 0; break;
case 1 : xi[0] = 1-a; xi[1] = a; break;
case 2 : xi[0] = 0; xi[1] = 1-a; break;
......@@ -525,7 +520,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
#else
double a = LQP[iqp][0];
double b = LQP[iqp][1];
switch(iebnd) {
switch(icl) {
case 0 : xi[0] = a; xi[1] = b; xi[2] = 0; break;
case 1 : xi[0] = 0; xi[1] = a; xi[2] = b; break;
case 2 : xi[0] = b; xi[1] = 0; xi[2] = a; break;
......@@ -546,11 +541,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
double c = 0;
double *vbnd = NULL;
if (n_value>0){
vbnd = malloc(sizeof(double)*n_value);
}
for (int i = 0; i < n_value; ++i){
for (int i = 0; i < wbnd->n_value; ++i){
vbnd[i] = 0;
}
double dc[D] = {0};
......@@ -561,8 +552,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
int node_local_id = bnd_node_local_id[el[i]];
if (node_local_id != -1)
for (int j = 0; j < n_value; j++){
vbnd[j] += v[node_local_id*n_value+j]*phi[i];
for (int j = 0; j < wbnd->n_value; j++){
vbnd[j] += v[node_local_id*wbnd->n_value+j]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
......@@ -580,7 +571,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
mu = problem->mu[0]*a + problem->mu[1]*(1-a);
}
const double jw = LQW[iqp]*detbnd;
f(problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,vbnd);
wbnd->cb(problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0 ? NULL :vbnd);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -589,12 +580,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
f(problem,n,g0,s,ds,c,dc,rho,mu,dt,iel,vbnd);
wbnd->cb(problem,n,g0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0? NULL :vbnd);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
f(problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,vbnd);
wbnd->cb(problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value ==0? NULL : vbnd);
ds[jfield*D+jd] = store;
}
int N2 = n_fields*N_SF;
......@@ -614,10 +605,9 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
}
free(vbnd);
}
}
free(x);
free(vbnd);
free(v);
}
free(bnd_node_local_id);
......@@ -786,14 +776,16 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys)
if (iphys == mesh->n_boundaries){
printf("Boundary tag \"%s\" not found.\n", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
exit(1);
}
for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){
forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
}
}
node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
......@@ -825,28 +817,32 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
free(all_local_vector);
}
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_fields, int nBnd, StrongBoundary *bnds)
static void apply_boundary_conditions(FluidProblem *problem)
{
for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
StrongBoundary *bnd = bnds + ibnd;
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys){
if (iphys == mesh->n_boundaries){
printf("Boundary tag \"%s\" not found.\n", bnd->tag);
exit(1);
}
double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
int j = mesh->phys_nodes[iphys][i];
MeshBoundary *mbnd = &problem->boundaries[iphys];
double *x = (double*)malloc(mbnd->n_nodes*sizeof(double)*D);
double *v = (double*)malloc(mbnd->n_nodes*sizeof(double));
for (int i = 0; i < mbnd->n_nodes; ++i){
int j = mbnd->nodes[i];
for (int k = 0; k < D; ++k)
x[i*D+k] = mesh->x[j*3+k];
}
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
bnd->apply(mbnd->n_nodes, x, v);
for (int i = 0; i < mbnd->n_nodes; ++i)
problem->solution[mbnd->nodes[i]*n_fields+bnd->field] = v[i];
free(x);
free(v);
}
......@@ -859,7 +855,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
struct timespec tic, toc;
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
apply_boundary_conditions(mesh, problem->solution, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
apply_boundary_conditions(problem);
double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
......@@ -997,6 +993,7 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
problem->particle_force = NULL;
problem->particle_volume = NULL;
problem->particle_velocity = NULL;
problem->boundaries = NULL;
return problem;
}
......@@ -1030,6 +1027,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->strong_boundaries);
mesh_tree_free(problem->mesh_tree);
mesh_free(problem->mesh);
mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
free(problem);
}
......@@ -1279,6 +1277,8 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
fluid_problem_compute_node_volume(problem);
if (problem->boundaries) mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
problem->boundaries = mesh_boundaries_create(problem->mesh);
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin)
......@@ -1477,7 +1477,7 @@ void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryC
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
void fluid_problem_set_open_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
......@@ -1489,11 +1489,7 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const cha
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
p->weak_boundaries[p->n_weak_boundaries-1].tag = strdup(tag);
f_cb *f = NULL;
if (strcasecmp(bndtype,"wall") == 0){
f = f_wall;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = 0;
}
else if (strcasecmp(bndtype,"velocity") == 0){
if (strcasecmp(bndtype,"velocity") == 0){
f = f_velocity;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = D+p->n_fluids-1;
}
......@@ -1529,22 +1525,16 @@ double *fluid_problem_old_porosity(const FluidProblem *p)
return p->oldporosity;
}
int fluid_problem_private_n_physical(const FluidProblem *p){
return p->mesh->n_phys;
int fluid_problem_private_n_boundaries(const FluidProblem *p){
return p->mesh->n_boundaries;
}
void fluid_problem_private_physical(const FluidProblem *p, int i, char **phys_name, int *phys_dim, int *phys_id, int *phys_n_nodes, int **phys_nodes, int *phys_n_boundaries, int **phys_boundaries)
void fluid_problem_boundary_name(const FluidProblem *p, int i, char **phys_name)
{
*phys_name = p->mesh->phys_name[i];
*phys_n_nodes = p->mesh->phys_n_nodes[i];
*phys_nodes = p->mesh->phys_nodes[i];
*phys_dim = p->mesh->phys_dimension[i];
*phys_id = p->mesh->phys_id[i];
*phys_n_boundaries = p->mesh->phys_n_boundaries[i];
*phys_boundaries = p->mesh->phys_boundaries[i];
*phys_name = p->mesh->boundary_names[i];
}
void fluid_problem_private_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements)
/*void fluid_problem_private_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements)
{
Mesh *m = malloc(sizeof(Mesh));
m->n_phys = 0;
......@@ -1562,9 +1552,9 @@ void fluid_problem_private_set_elements(FluidProblem *p, int n_nodes, double *x,
m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_elements);
memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);
fluid_problem_set_mesh(p,m);
}
}*/
void fluid_problem_private_add_physical(FluidProblem *p, int tag, int dim, const char *name, int n_nodes, int *nodes, int n_bnd, int *bnd)
/*void fluid_problem_private_add_physical(FluidProblem *p, int tag, int dim, const char *name, int n_nodes, int *nodes, int n_bnd, int *bnd)
{
Mesh *m = p->mesh;
int i = m->n_phys;
......@@ -1585,5 +1575,5 @@ void fluid_problem_private_add_physical(FluidProblem *p, int tag, int dim, const
m->phys_boundaries = realloc(m->phys_boundaries,m->n_phys*sizeof(int*));
m->phys_boundaries[i] = malloc(n_bnd*sizeof(int)*2);
memcpy(m->phys_boundaries[i],bnd,n_bnd*sizeof(int)*2);
}
}*/
......@@ -59,6 +59,7 @@ struct FluidProblem {
double coeffStab;
Mesh *mesh;
MeshTree *mesh_tree;
MeshBoundary *boundaries;
double *porosity;
double *oldporosity;
double *concentration;
......@@ -103,7 +104,7 @@ int fluid_problem_n_fields(const FluidProblem *p);
double *fluid_problem_solution(const FluidProblem *p);
double *fluid_problem_coordinates(const FluidProblem *p);
int fluid_problem_n_nodes(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply);
void fluid_problem_set_open_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply);
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
......
......@@ -27,6 +27,17 @@
#include <string.h>
#include "tools.h"
#if DIMENSION == 2
static int elbnd[][2] = {{0,1},{1,2},{2,0},{1,0},{2,1},{0,2}};
#else
static int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
typedef struct {
int tag;
int node[DIMENSION];
} BoundaryElement;
int edgecmp(const void *pe0, const void *pe1) {
int *e0 = (int*)pe0;
int *e1 = (int*)pe1;
......@@ -34,26 +45,49 @@ int edgecmp(const void *pe0, const void *pe1) {
return e0[1]-e1[1];
}
static void mesh_gen_edges(Mesh *m) {
static void swapint(int *a, int *b) {
int c = *a;
*a = *b;
*b = c;
}
static void mesh_gen_edges(Mesh *m, int n_bnd_elements, BoundaryElement *bnd_elements) {
#if DIMENSION == 2
int n_half_edges = m->n_elements*3;
int *halfedges = malloc(sizeof(int)*n_half_edges*4);
int n_half_edges = m->n_elements*3+n_bnd_elements;
int *halfedges = malloc(sizeof(int)*n_half_edges*4); //n0 n1 {(el,closure),(-1,tag)}
int edgemap[3][2] = {{0,1},{1,2},{2,0}};
int id = 0;
for (int i=0; i<m->n_elements; ++i) {
for (int j = 0; j < 3; ++j) {
int n0 = m->elements[i*3+edgemap[j][0]], n1 = m->elements[i*3+edgemap[j][1]];
int n0 = m->elements[i*3+edgemap[j][0]];
int n1 = m->elements[i*3+edgemap[j][1]];
if (n0 < n1) {
halfedges[(i*3+j)*4+0] = n0;
halfedges[(i*3+j)*4+1] = n1;
halfedges[(i*3+j)*4+3] = j;
halfedges[id*4+0] = n0;
halfedges[id*4+1] = n1;
halfedges[id*4+3] = j;
}
else {
halfedges[(i*3+j)*4+0] = n1;
halfedges[(i*3+j)*4+1] = n0;
halfedges[(i*3+j)*4+3] = j+3;
halfedges[id*4+0] = n1;
halfedges[id*4+1] = n0;
halfedges[id*4+3] = j+3;
}
halfedges[(i*3+j)*4+2] = i;
halfedges[id*4+2] = i;
id++;
}
}
for (int i=0; i<n_bnd_elements; ++i) {
BoundaryElement *e = bnd_elements + i;
if (e->node[0] < e->node[1]) {
halfedges[id*4+0] = e->node[0];
halfedges[id*4+1] = e->node[1];
}
else {
halfedges[id*4+0] = e->node[1];
halfedges[id*4+1] = e->node[0];
}
halfedges[id*4+2] = -1;
halfedges[id*4+3] = e->tag;
id ++;
}
qsort(halfedges,n_half_edges,4*sizeof(int),edgecmp);
int n_edges = 0;
......@@ -79,13 +113,19 @@ static void mesh_gen_edges(Mesh *m) {
m->interfaces = edges;
free(halfedges);
for (int i = 0; i < n_edges; ++i) {
if (edges[i*4+1] >=3)
// ensure boundary tag is on the second part of the closure
if (edges[i*4+0] < 0) {
swapint(edges+i*4,edges+i*4+2);
swapint(edges+i*4+1,edges+i*4+3);
}
// ensure edge is alligned with first element
if (edges[i*4+1] >= 3)
edges[i*4+1] -= 3;
if (edges[i*4+3] >= 0 && edges[i*4+3] < 3)
if (edges[i*4+2] >= 0 && edges[i*4+3] < 3)
edges[i*4+3] += 3;
}
#else
// todo Matthieu
// todo Jon
#endif
}
......@@ -122,22 +162,6 @@ static int read_quoted_string(FILE *f, char *w, size_t maxl) {
return 0;
}
static void mesh_add_physical_node(Mesh *m, int dim, int id, int node)
{
int i = 0;
for (i = 0; i < m->n_phys; ++i) {
if(m->phys_id[i] == id && m->phys_dimension[i] == dim) {
if((m->phys_n_nodes[i]/2)*2 == m->phys_n_nodes[i])
m->phys_nodes[i] = realloc(m->phys_nodes[i], 2*(m->phys_n_nodes[i]+1)*sizeof(int));
m->phys_nodes[i][m->phys_n_nodes[i]] = node;
m->phys_n_nodes[i] += 1;
break;
}
}
if (i == m->n_phys)
printf("Unknown physical tag %i with dimension %i.\n", id, dim);
}
static int int_cmp(const void *p0, const void *p1)
{
int i0 = *(int*)p0;
......@@ -173,23 +197,6 @@ static void sortDi(int i[DIMENSION]) {
#endif
}
typedef struct {
int tag;
int node[DIMENSION];
} BoundaryElement;
int bnd_element_cmp(const void *ps0, const void *ps1) {
BoundaryElement *s0 = (BoundaryElement*) ps0;
BoundaryElement *s1 = (BoundaryElement*) ps1;
if (s0->node[0] != s1->node[0]) return s0->node[0]-s1->node[0];
#if DIMENSION == 2
return s0->node[1]-s1->node[1];
#else
if (s0->node[1] != s1->node[1]) return s0->node[1]-s1->node[1];
return s0->node[2]-s1->node[2];
#endif
}
Mesh *mesh_load_msh(const char *filename)
{
Mesh *m = malloc(sizeof(Mesh));
......@@ -212,26 +219,29 @@ Mesh *mesh_load_msh(const char *filename)
}
check_word_in_file(f,"$EndMeshFormat");
check_word_in_file(f,"$PhysicalNames");
if (fscanf(f, "%d", &m->n_phys) != 1){
int nphysfile;
if (fscanf(f, "%d", &nphysfile) != 1){
printf("Cannot read physical names\n");
return NULL;
}
m->phys_name = malloc(sizeof(char*)*m->n_phys);
m->phys_nodes = malloc(sizeof(int*)*m->n_phys);
m->phys_dimension = malloc(sizeof(int)*m->n_phys);
m->phys_n_nodes = malloc(sizeof(int)*m->n_phys);
m->phys_id = malloc(sizeof(int)*m->n_phys);
for (int i = 0; i < m->n_phys; ++i) {
if (fscanf(f, "%i %i", &m->phys_dimension[i], &m->phys_id[i]) != 2){
m->boundary_names = malloc(sizeof(char*)*nphysfile);
int *physid = malloc(sizeof(int)*nphysfile);
int nphys = 0;
for (int i = 0; i < nphysfile; ++i) {
int physdim, pid;
if (fscanf(f, "%i %i", &physdim, &pid) != 2){
printf("Cannot read physical names\n");
return NULL;
}
char pname[256];
if(read_quoted_string(f, pname, 256) < 0) return NULL;
m->phys_name[i] = strdup(pname);
m->phys_nodes[i] = NULL;
m->phys_n_nodes[i] = 0;
if (physdim == DIMENSION-1) {
physid[nphys] = pid;
m->boundary_names[nphys] = strdup(pname);
nphys++;
}
}
m->n_boundaries = nphys;
check_word_in_file(f,"$EndPhysicalNames");
check_word_in_file(f,"$Nodes");
if (fscanf(f, "%d", &m->n_nodes) != 1){
......@@ -264,7 +274,7 @@ Mesh *mesh_load_msh(const char *filename)
m->n_elements = 0;
for (int i = 0; i < n_el; ++i) {
int eid, etype, nflags, flag;
int phys=0;
int physfile=0;
if (fscanf(f, "%d %d %d", &eid, &etype, &nflags) != 3){
printf("Cannot read elements.\n");
return NULL;
......@@ -275,9 +285,10 @@ Mesh *mesh_load_msh(const char *filename)
return NULL;
}
if (j == 0)
phys = flag;
physfile = flag;
}
int n[DIMENSION+1];
int phys = -1;
switch (etype) {
case(1):
if (fscanf(f, "%d %d", &n[0], &n[1]) != 2){
......@@ -285,20 +296,23 @@ Mesh *mesh_load_msh(const char *filename)
return NULL;
}
#if DIMENSION==2
bnd_elements[n_bnd_elements].tag = phys;
for (int iD = 0; iD < 2; ++iD) bnd_elements[n_bnd_elements].node[iD] = n[iD]-1;
sortDi(bnd_elements[n_bnd_elements].node);
n_bnd_elements++;
phys = -1;
for (int in = 0; in < m->n_boundaries; ++in) {
if (physid[in] == physfile) phys = in;
}
if (phys != -1) {
bnd_elements[n_bnd_elements].tag = phys;
for (int iD = 0; iD < 2; ++iD) bnd_elements[n_bnd_elements].node[iD] = n[iD]-1;
sortDi(bnd_elements[n_bnd_elements].node);
n_bnd_elements++;
}
#endif
mesh_add_physical_node(m, 1, phys, n[0]-1);
mesh_add_physical_node(m, 1, phys, n[1]-1);
break;
case(15):
if (fscanf(f, "%d", &n[0]) != 1){