diff --git a/opm/core/pressure/flow_bc.c b/opm/core/pressure/flow_bc.c index dba5701ff..82451f947 100644 --- a/opm/core/pressure/flow_bc.c +++ b/opm/core/pressure/flow_bc.c @@ -19,6 +19,7 @@ #include #include +#include #include @@ -50,48 +51,63 @@ static void initialise_structure(struct FlowBoundaryConditions *fbc) /* ---------------------------------------------------------------------- */ { - fbc->nbc = 0; - fbc->cpty = 0; + fbc->nbc = 0; + fbc->cond_cpty = 0; + fbc->face_cpty = 0; - fbc->type = NULL; - fbc->value = NULL; - fbc->face = NULL; + fbc->cond_pos = NULL; + fbc->type = NULL; + fbc->value = NULL; + fbc->face = NULL; } /* ---------------------------------------------------------------------- */ static int expand_tables(size_t nbc, + size_t nf , struct FlowBoundaryConditions *fbc) /* ---------------------------------------------------------------------- */ { - int ok; + int ok_cond, ok_face; size_t alloc_sz; - void *p1, *p2, *p3; + void *p1, *p2, *p3, *p4; - ok = nbc <= fbc->cpty; + ok_cond = nbc <= fbc->cond_cpty; + ok_face = nf <= fbc->face_cpty; - if (! ok) { - alloc_sz = alloc_size(nbc, fbc->cpty); + if (! ok_cond) { + alloc_sz = alloc_size(nbc, fbc->cond_cpty); - p1 = realloc(fbc->type , alloc_sz * sizeof *fbc->type ); - p2 = realloc(fbc->value, alloc_sz * sizeof *fbc->value); - p3 = realloc(fbc->face , alloc_sz * sizeof *fbc->face ); + p1 = realloc(fbc->type , (alloc_sz + 0) * sizeof *fbc->type ); + p2 = realloc(fbc->value , (alloc_sz + 0) * sizeof *fbc->value ); + p3 = realloc(fbc->cond_pos, (alloc_sz + 1) * sizeof *fbc->cond_pos); - ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL); + ok_cond = (p1 != NULL) && (p2 != NULL) && (p3 != NULL); - if (ok) { - fbc->type = p1; - fbc->value = p2; - fbc->face = p3; + if (p1 != NULL) { fbc->type = p1; } + if (p2 != NULL) { fbc->value = p2; } + if (p3 != NULL) { fbc->cond_pos = p3; } - fbc->cpty = alloc_sz; - } else { - free(p3); free(p2); free(p1); + if (ok_cond) { + fbc->cond_cpty = alloc_sz; } } - return ok; + if (! ok_face) { + alloc_sz = alloc_size(nf, fbc->face_cpty); + + p4 = realloc(fbc->face, alloc_sz * sizeof *fbc->face); + + ok_face = p4 != NULL; + + if (ok_face) { + fbc->face = p4; + fbc->face_cpty = alloc_sz; + } + } + + return ok_cond && ok_face; } @@ -116,12 +132,14 @@ flow_conditions_construct(size_t nbc) if (fbc != NULL) { initialise_structure(fbc); - ok = expand_tables(nbc, fbc); + ok = expand_tables(nbc, nbc, fbc); if (! ok) { flow_conditions_destroy(fbc); fbc = NULL; + } else { + fbc->cond_pos[0] = 0; } } @@ -138,9 +156,10 @@ flow_conditions_destroy(struct FlowBoundaryConditions *fbc) /* ---------------------------------------------------------------------- */ { if (fbc != NULL) { - free(fbc->face ); - free(fbc->value); - free(fbc->type ); + free(fbc->face ); + free(fbc->cond_pos); + free(fbc->value ); + free(fbc->type ); } free(fbc); @@ -159,14 +178,38 @@ flow_conditions_append(enum FlowBCType type , struct FlowBoundaryConditions *fbc ) /* ---------------------------------------------------------------------- */ { - int ok; + return flow_conditions_append_multi(type, 1, &face, value, fbc); +} - ok = expand_tables(fbc->nbc + 1, fbc); + +/* ---------------------------------------------------------------------- */ +/* Append a new boundary condition that affects multiple interfaces. + * + * Return one (1) if successful, and zero (0) otherwise. */ +/* ---------------------------------------------------------------------- */ +int +flow_conditions_append_multi(enum FlowBCType type , + size_t nfaces, + const int *faces , + double value , + struct FlowBoundaryConditions *fbc ) +/* ---------------------------------------------------------------------- */ +{ + int ok; + size_t nbc; + + nbc = fbc->nbc; + + ok = expand_tables(nbc + 1, fbc->cond_pos[ nbc ] + nfaces, fbc); if (ok) { - fbc->type [ fbc->nbc ] = type ; - fbc->value[ fbc->nbc ] = value; - fbc->face [ fbc->nbc ] = face ; + memcpy(fbc->face + fbc->cond_pos[ nbc ], + faces, nfaces * sizeof *faces); + + fbc->type [ nbc ] = type; + fbc->value[ nbc ] = value; + + fbc->cond_pos[ nbc + 1 ] = fbc->cond_pos[ nbc ] + nfaces; fbc->nbc += 1; } diff --git a/opm/core/pressure/flow_bc.h b/opm/core/pressure/flow_bc.h index 7d209db81..65466d3e1 100644 --- a/opm/core/pressure/flow_bc.h +++ b/opm/core/pressure/flow_bc.h @@ -38,7 +38,11 @@ enum FlowBCType { BC_NOFLOW , * The field 'cpty' is for internal use by the implementation. */ struct FlowBoundaryConditions { size_t nbc; /* Current number of bdry. conditions */ - size_t cpty; /* Internal management. Do not touch */ + + size_t face_cpty; /* Internal management. Do not touch */ + size_t cond_cpty; /* Internal management. Do not touch */ + + size_t *cond_pos; /* Indirection pointer into '.face' */ enum FlowBCType *type; /* Condition type */ double *value; /* Condition value (target) */ @@ -67,6 +71,16 @@ flow_conditions_append(enum FlowBCType type , double value, struct FlowBoundaryConditions *fbc ); +/* Append a new boundary condition that affects multiple interfaces. + * + * Return one (1) if successful, and zero (0) otherwise. */ +int +flow_conditions_append_multi(enum FlowBCType type , + size_t nfaces, + const int *faces , + double value , + struct FlowBoundaryConditions *fbc ); + /* Clear existing set of boundary conditions */ void diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index 56bcb995b..b40bde83c 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -161,7 +161,7 @@ assemble_bc_contrib(struct UnstructuredGrid *G , int is_neumann; int i, f, c1, c2; - size_t j; + size_t j, ix; is_neumann = 1; @@ -169,17 +169,19 @@ assemble_bc_contrib(struct UnstructuredGrid *G , if (bc->type[ i ] == BC_PRESSURE) { is_neumann = 0; - f = bc->face[ i ]; - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; + for (j = bc->cond_pos[ i ]; j < bc->cond_pos[i + 1]; j++) { + f = bc->face[ j ]; + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; - assert ((c1 < 0) ^ (c2 < 0)); /* BCs on ext. faces only */ + assert ((c1 < 0) ^ (c2 < 0)); /* BCs on ext. faces only */ - c1 = (c1 >= 0) ? c1 : c2; - j = csrmatrix_elm_index(c1, c1, h->A); + c1 = (c1 >= 0) ? c1 : c2; + ix = csrmatrix_elm_index(c1, c1, h->A); - h->A->sa[ j ] += trans[ f ]; - h->b [ c1] += trans[ f ] * bc->value[ i ]; + h->A->sa[ ix ] += trans[ f ]; + h->b [ c1 ] += trans[ f ] * bc->value[ i ]; + } } /* Other types currently not handled */