自学教程：C++ ERR_CheckGo函数代码示例

int32 conn_alloc(MeshConnectivity *conn, uint32 num, uint32 n_incident){  int32 ret = RET_OK;  if ((conn->num > 0) && (conn->n_incident > 0)) {    conn_free(conn);  }  if (num > 0) {    conn->num = num;    conn->offsets = alloc_mem(uint32, num + 1);    ERR_CheckGo(ret);  }  if (n_incident > 0) {    conn->n_incident = n_incident;    conn->indices = alloc_mem(uint32, n_incident);    ERR_CheckGo(ret);  } else if (num == 0) { // Empty connectivity.    conn->n_incident = n_incident;    conn->indices = 0;  } end_label:  if (ERR_Chk) {    conn_free(conn);  }  return(ret);}

int32 mem_check_ptr(char *p, int lineNo, char *funName,                    char *fileName, char *dirName){  int32 ret = RET_OK;  float64 *endptr;  size_t hsize = sizeof(AllocSpaceAlign);  AllocSpace *head = (AllocSpace *) (p - hsize);  if (head->cookie != AL_CookieValue) {    errput("%s, %s, %s, %d: ptr: %p, cookie: %d/n",           dirName, fileName, funName, lineNo,           p, head->cookie);    if (head->cookie == AL_AlreadyFreed) {      errput("memory was already freed!/n");    }    ERR_CheckGo(ret);  }  endptr = (float64 *) (p + head->size);  if (endptr[0] != AL_CookieValue) {    errput("%s %s %s %d:/n",           dirName, fileName, funName, lineNo);    if (endptr[0] == AL_AlreadyFreed) {      errput("already freed!/n");    } else {      errput("damaged tail!/n");    }    ERR_CheckGo(ret);  } end_label:  return(ret);}

/*!  @par Revision history:  - 01.02.2008, c*/int32 dw_mass_scalar( FMField *out, FMField *coef,		      FMField *state, FMField *bf, VolumeGeometry *vg,		      int32 isDiff ){  int32 ii, dim, nQP, nEP, ret = RET_OK;  FMField *ftfp = 0, *ftf = 0, *cftf = 0;  nQP = vg->bfGM->nLev;  nEP = vg->bfGM->nCol;  dim = vg->bfGM->nRow;  if (isDiff) {    fmf_createAlloc( &ftf, 1, nQP, nEP, nEP );    fmf_createAlloc( &cftf, 1, nQP, nEP, nEP );    fmf_mulATB_nn( ftf, bf, bf );    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( vg->det, ii );      if (coef->nCell > 1) {	FMF_SetCell( coef, ii );      }      fmf_mulAF( cftf, ftf, coef->val );      fmf_sumLevelsMulF( out, cftf, vg->det->val );      ERR_CheckGo( ret );    }  } else {    fmf_createAlloc( &ftfp, 1, nQP, nEP, 1 );    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( state, ii );      FMF_SetCell( vg->det, ii );      if (coef->nCell > 1) {	FMF_SetCell( coef, ii );      }      bf_actt( ftfp, bf, state );      fmf_mul( ftfp, coef->val );      fmf_sumLevelsMulF( out, ftfp, vg->det->val );      ERR_CheckGo( ret );    }  } end_label:  if (isDiff) {    fmf_freeDestroy( &ftf );    fmf_freeDestroy( &cftf );  } else {    fmf_freeDestroy( &ftfp );  }  return( ret );}

/*!  @par Revision history:  - 09.03.2009, c*/int32 dw_surf_mass_scalar( FMField *out, FMField *coef,			   FMField *state, FMField *bf, SurfaceGeometry *sg,			   int32 isDiff ){  int32 ii, nFP, ret = RET_OK;  FMField *ftfp = 0, *ftf = 0, *cftf = 0;  nFP = bf->nCol;  if (isDiff) {    fmf_createAlloc( &ftf, 1, sg->nQP, nFP, nFP );    fmf_createAlloc( &cftf, 1, sg->nQP, nFP, nFP );    fmf_mulATB_nn( ftf, bf, bf );    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( sg->det, ii );      if (coef->nCell > 1) {	FMF_SetCell( coef, ii );      }      fmf_mulAF( cftf, ftf, coef->val );      fmf_sumLevelsMulF( out, cftf, sg->det->val );      ERR_CheckGo( ret );    }  } else {    fmf_createAlloc( &ftfp, 1, sg->nQP, nFP, 1 );    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( sg->det, ii );      if (coef->nCell > 1) {	FMF_SetCell( coef, ii );      }      FMF_SetCell( state, ii );      bf_actt( ftfp, bf, state );      fmf_mulAF( ftfp, ftfp, coef->val );      fmf_sumLevelsMulF( out, ftfp, sg->det->val );      ERR_CheckGo( ret );    }  }end_label:  if (isDiff) {    fmf_freeDestroy( &ftf );    fmf_freeDestroy( &cftf );  } else {    fmf_freeDestroy( &ftfp );  }  return( ret );}

int32 mesh_transpose(Mesh *mesh, int32 d1, int32 d2){  int32 ret = RET_OK;  uint32 n_incident;  uint32 ii;  uint32 *nd2 = 0;  uint32 D = mesh->topology->max_dim;  MeshEntityIterator it2[1], it1[1];  MeshConnectivity *c12 = 0; // d1 -> d2 - to compute  debprintf("transpose %d -> %d/n", d1, d2);  if (d1 >= d2) {    errput("d1 must be smaller than d2 in mesh_transpose()!/n");    ERR_CheckGo(ret);  }  c12 = mesh->topology->conn[IJ(D, d1, d2)];  // Count entities of d2 -> d1.  conn_alloc(c12, mesh->topology->num[d1], 0);  ERR_CheckGo(ret);  nd2 = c12->offsets + 1;  for (mei_init(it2, mesh, d2); mei_go(it2); mei_next(it2)) {    for (mei_init_conn(it1, it2->entity, d1); mei_go(it1); mei_next(it1)) {      nd2[it1->entity->ii]++;    }  }  // c12->offsets now contains counts - make a cumsum to get offsets.  for (ii = 1; ii < c12->num + 1; ii++) {    c12->offsets[ii] += c12->offsets[ii-1];  }  n_incident = c12->offsets[c12->num];  debprintf("transpose n_incident (%d -> %d): %d/n", d1, d2, n_incident);  // Fill in the indices.  conn_alloc(c12, 0, n_incident);  ERR_CheckGo(ret);  for (ii = 0; ii < c12->n_incident; ii++) {    c12->indices[ii] = UINT32_None; // "not set" value.  }  for (mei_init(it2, mesh, d2); mei_go(it2); mei_next(it2)) {    for (mei_init_conn(it1, it2->entity, d1); mei_go(it1); mei_next(it1)) {      conn_set_to_free(c12, it1->entity->ii, it2->entity->ii);      ERR_CheckGo(ret);    }  } end_label:  return(ret);}

/*!  @par Revision history:  - 23.03.2007, c*/int32 d_of_nsSurfMinDPress( FMField *out, FMField *pressure,                            float64 weight, float64 bpress,			    SurfaceGeometry *sg, int32 isDiff ){  int32 ii, iqp, nQP, ret = RET_OK;  float64 aux;  FMField *out1 = 0, *pressQP = 0;  nQP = sg->det->nLev;  if (isDiff == 0) {    fmf_createAlloc( &out1, 1, 1, 1, 1 );    fmf_createAlloc( &pressQP, 1, nQP, 1, 1 );    aux = 0.0;    for (ii = 0; ii < pressure->nCell; ii++) {      FMF_SetCell( pressure, ii );      FMF_SetCell( sg->det, ii );      for (iqp = 0; iqp < nQP; iqp++) {	pressQP->val[iqp] -= pressure->val[iqp] - bpress;      }      fmf_sumLevelsMulF( out1, pressQP, sg->det->val );      aux += out1->val[0];      ERR_CheckGo( ret );    }    out->val[0] = aux * weight;  } else {    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( sg->det, ii );      FMF_SetCellX1( sg->bf, ii );      fmf_sumLevelsTMulF( out, sg->bf, sg->det->val );      ERR_CheckGo( ret );    }    fmfc_mulC( out, weight );  } end_label:  if (isDiff == 0) {    fmf_freeDestroy( &out1 );    fmf_freeDestroy( &pressQP );  }  return( ret );}

/*!  @par Revision history:  - 23.01.2006, c  - 02.03.2006  - 27.07.2006*/int32 d_of_nsMinGrad( FMField *out, FMField *grad,		      FMField *viscosity, VolumeGeometry *vg ){  int32 ii, nQP, ret = RET_OK;  float64 aux;  FMField *out1 = 0, *gvel2 = 0;  nQP = vg->bfGM->nLev;  fmf_createAlloc( &out1, 1, 1, 1, 1 );  fmf_createAlloc( &gvel2, 1, nQP, 1, 1 );  FMF_SetFirst( out );  aux = 0.0;  for (ii = 0; ii < grad->nCell; ii++) {    FMF_SetCell( grad, ii );    FMF_SetCell( viscosity, ii );    FMF_SetCell( vg->det, ii );    fmf_mulATB_nn( gvel2, grad, grad );    fmf_mul( gvel2, viscosity->val );    fmf_sumLevelsMulF( out1, gvel2, vg->det->val );    aux += out1->val[0];    ERR_CheckGo( ret );  }  out->val[0] = aux * 0.5; end_label:  fmf_freeDestroy( &out1 );  fmf_freeDestroy( &gvel2 );  return( ret );}

/*!  @par Revision history:  - c: 25.03.2008*/int32 de_cauchy_stress( FMField *out, FMField *strain,			FMField *mtxD,  VolumeGeometry *vg,			int32 mode ){  int32 ii, dim, sym, nQP, ret = RET_OK;  FMField *stress = 0;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  sym = (dim + 1) * dim / 2;  fmf_createAlloc( &stress, 1, nQP, sym, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( mtxD, ii );    FMF_SetCell( strain, ii );    FMF_SetCell( vg->det, ii );    fmf_mulAB_nn( stress, mtxD, strain );    fmf_sumLevelsMulF( out, stress, vg->det->val );    if (mode == 1) {      FMF_SetCell( vg->volume, ii );      fmf_mulC( out, 1.0 / vg->volume->val[0] );    }    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &stress );  return( ret );}

/*!  @par Revision history:  - 21.09.2006, c*/int32 de_cauchy_strain( FMField *out, FMField *strain,			VolumeGeometry *vg, int32 mode ){  int32 ii, dim, sym, nQP, ret = RET_OK;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  sym = (dim + 1) * dim / 2;  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( strain, ii );    FMF_SetCell( vg->det, ii );    fmf_sumLevelsMulF( out, strain, vg->det->val );    if (mode == 1) {      FMF_SetCell( vg->volume, ii );      fmf_mulC( out, 1.0 / vg->volume->val[0] );    }    ERR_CheckGo( ret );  } end_label:  return( ret );}

int32 dw_lin_prestress( FMField *out, FMField *stress, VolumeGeometry *vg,			int32 *elList, int32 elList_nRow, int32 isDiff ){  int32 ii, iel, dim, nQP, nEP, ret = RET_OK;  FMField *res = 0;  nQP = vg->bfGM->nLev;  nEP = vg->bfGM->nCol;  dim = vg->bfGM->nRow;  fmf_createAlloc( &res, 1, nQP, dim * nEP, 1 );    for (ii = 0; ii < elList_nRow; ii++) {      iel = elList[ii];      FMF_SetCell( out, ii );      FMF_SetCell( vg->bfGM, iel );      FMF_SetCell( vg->det, iel );      FMF_SetCell( stress, iel );      form_sdcc_actOpGT_VS3( res, vg->bfGM, stress );      fmf_sumLevelsMulF( out, res, vg->det->val );      ERR_CheckGo( ret );    } end_label:    fmf_freeDestroy( &res );   return( ret );}

int32 d_surf_lcouple( FMField *out, FMField *stateP, FMField *gradQ, FMField *coef,                      Mapping *sg ){    int32 ii, dim, nQP, ret = RET_OK;    FMField *aux1 = 0, *aux2 = 0;    nQP = stateP->nLev;    dim = coef->nRow;    fmf_createAlloc(&aux1, 1, nQP, dim, 1);    fmf_createAlloc(&aux2, 1, nQP, 1, 1);    for (ii = 0; ii < out->nCell; ii++) {        if (coef->nCell > 1)            FMF_SetCell( coef, ii );        FMF_SetCell( out, ii );        FMF_SetCell( sg->det, ii );        FMF_SetCell(stateP, ii);        FMF_SetCell(gradQ, ii);        fmf_mulAB_nn(aux1, coef, stateP);        fmf_mulATB_nn(aux2, gradQ, aux1);        fmf_sumLevelsMulF(out, aux2, sg->det->val);        ERR_CheckGo( ret );    }end_label:    fmf_freeDestroy(&aux1);    fmf_freeDestroy(&aux2);    return( ret );}

int32 dw_lin_strain_fib( FMField *out, FMField *mtxD, FMField *mat,			 Mapping *vg ){  int32 ii, dim, sym, nQP, nEP, ret = RET_OK;  FMField *aux1 = 0, *aux2 = 0;  nQP = vg->bfGM->nLev;  nEP = vg->bfGM->nCol;  dim = vg->bfGM->nRow;  sym = (dim + 1) * dim / 2;  fmf_createAlloc( &aux1, 1, nQP, nEP * dim, sym );  fmf_createAlloc( &aux2, 1, nQP, nEP * dim, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( mtxD, ii );    FMF_SetCell( mat, ii );    FMF_SetCell( vg->bfGM, ii );    FMF_SetCell( vg->det, ii );    form_sdcc_actOpGT_M3( aux1, vg->bfGM, mtxD );    fmf_mulAB_nn( aux2, aux1, mat );    fmf_sumLevelsMulF( out, aux2, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &aux1 );  fmf_freeDestroy( &aux2 );  return( ret );}

int32 dw_lin_prestress( FMField *out, FMField *stress, Mapping *vg ){  int32 ii, dim, nQP, nEP, ret = RET_OK;  FMField *res = 0;  nQP = vg->bfGM->nLev;  nEP = vg->bfGM->nCol;  dim = vg->bfGM->nRow;  fmf_createAlloc( &res, 1, nQP, dim * nEP, 1 );    for (ii = 0; ii < out->nCell; ii++) {      FMF_SetCell( out, ii );      FMF_SetCell( vg->bfGM, ii );      FMF_SetCell( vg->det, ii );      FMF_SetCell( stress, ii );      form_sdcc_actOpGT_VS3( res, vg->bfGM, stress );      fmf_sumLevelsMulF( out, res, vg->det->val );      ERR_CheckGo( ret );    } end_label:    fmf_freeDestroy( &res );  return( ret );}

/*!  mode == 0: /int pq  mode == 1: /int pq div V  Works for both scalars and vectors.  @par Revision history:  - 24.02.2006, c  - 27.02.2006*/int32 d_sd_volume_dot( FMField *out, FMField *stateP, FMField *stateQ,                       FMField *divMV, VolumeGeometry *vg, int32 mode ){  int32 ii, nQP, ret = RET_OK;  FMField *pq = 0;  nQP = vg->bfGM->nLev;  fmf_createAlloc( &pq, 1, nQP, 1, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( stateP, ii );    FMF_SetCell( stateQ, ii );    FMF_SetCell( vg->det, ii );    fmf_mulATB_nn( pq, stateP, stateQ );    if (mode == 1) {      FMF_SetCell( divMV, ii );      fmf_mul( pq, divMV->val );    }    fmf_sumLevelsMulF( out, pq, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &pq );  return( ret );}

/*!  @par Revision history:  - 22.03.2006, c*/int32 d_sd_convect( FMField *out, FMField *stateU, FMField *gradU,		    FMField *stateW, FMField *divMV, FMField *gradMV,		    VolumeGeometry *vg_u, int32 mode ){  int32 ii, dim, nQP, ret = RET_OK;  FMField *aux11 = 0, *aux = 0, *fwgu = 0, *fwgugmv = 0;  FMField gum[1], gmvm[1];  nQP = vg_u->bfGM->nLev;  dim = vg_u->bfGM->nRow;  gum->nAlloc = -1;  fmf_pretend( gum, gradU->nCell, nQP, dim, dim, gradU->val0 );  fmf_createAlloc( &fwgu, 1, nQP, 1, dim );  fmf_createAlloc( &aux11, 1, nQP, 1, 1 );  if (mode == 1) {    gmvm->nAlloc = -1;    fmf_pretend( gmvm, gradMV->nCell, nQP, dim, dim, gradMV->val0 );    fmf_createAlloc( &fwgugmv, 1, nQP, 1, dim );    fmf_createAlloc( &aux, 1, nQP, 1, 1 );  }  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( stateU, ii );    FMF_SetCell( stateW, ii );    FMF_SetCell( gum, ii );    FMF_SetCell( vg_u->det, ii );    fmf_mulATB_nn( fwgu, stateW, gum );    fmf_mulAB_nn( aux11, fwgu, stateU );    if (mode == 1) {      FMF_SetCell( divMV, ii );      FMF_SetCell( gmvm, ii );      fmf_mul( aux11, divMV->val );      fmf_mulAB_nn( fwgugmv, fwgu, gmvm );      fmf_mulAB_nn( aux, fwgugmv, stateU );      fmf_subAB_nn( aux11, aux11, aux );    }    fmf_sumLevelsMulF( out, aux11, vg_u->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &aux11 );  fmf_freeDestroy( &fwgu );  if (mode == 1) {    fmf_freeDestroy( &fwgugmv );    fmf_freeDestroy( &aux );  }  return( ret );}

/*!  @par Revision history:  - 04.09.2007, c*/int32 d_mass_scalar( FMField *out, FMField *coef,		     FMField *stateP, FMField *stateQ,		     FMField *bf, VolumeGeometry *vg ){  int32 ii, dim, nQP, ret = RET_OK;  FMField *qftfp = 0;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  fmf_createAlloc( &qftfp, 1, nQP, 1, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( stateP, ii );    FMF_SetCell( stateQ, ii );    FMF_SetCell( vg->det, ii );    if (coef->nCell > 1) {      FMF_SetCell( coef, ii );    }    fmf_mulATB_nn( qftfp, stateQ, stateP );    fmf_mul( qftfp, coef->val );    fmf_sumLevelsMulF( out, qftfp, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &qftfp );  return( ret );}

/*!  @par Revision history:  - c: 05.03.2008, r: 05.03.2008*/int32 d_biot_div( FMField *out, float64 coef, FMField *state, FMField *strain,		  FMField *mtxD, Mapping *vg ){  int32 ii, nQP, ret = RET_OK;  FMField *dtgu = 0, *pftdtgu = 0;  nQP = vg->bfGM->nLev;/*   fmf_print( mtxD, stdout, 0 ); */  fmf_createAlloc( &dtgu, 1, nQP, 1, 1 );  fmf_createAlloc( &pftdtgu, 1, nQP, 1, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( mtxD, ii );    FMF_SetCell( vg->det, ii );    FMF_SetCell( state, ii );    FMF_SetCell( strain, ii );    fmf_mulATB_nn( dtgu, mtxD, strain );    fmf_mulATB_nn( pftdtgu, state, dtgu );    fmf_sumLevelsMulF( out, pftdtgu, vg->det->val );    ERR_CheckGo( ret );  }  // E.g. 1/dt.  fmfc_mulC( out, coef ); end_label:  fmf_freeDestroy( &dtgu );  fmf_freeDestroy( &pftdtgu );  return( ret );}

/*!  @par Revision history:  - 30.07.2007, from dw_hdpm_cache()*/int32 dq_cauchy_strain( FMField *out, FMField *state, int32 offset,			VolumeGeometry *vg,			int32 *conn, int32 nEl, int32 nEP ){  int32 ii, dim, sym, nQP, ret = RET_OK;  FMField *st = 0, *disG = 0;  state->val = FMF_PtrFirst( state ) + offset;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  sym = (dim + 1) * dim / 2;  fmf_createAlloc( &st, 1, 1, nEP, dim );  fmf_createAlloc( &disG, 1, nQP, dim, dim );  for (ii = 0; ii < nEl; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( vg->bfGM, ii );    ele_extractNodalValuesNBN( st, state, conn + nEP * ii );    fmf_mulAB_n1( disG, vg->bfGM, st );    form_sdcc_strainCauchy_VS( out, disG );/*       fmf_print( out, stdout, 0 ); *//*       sys_pause(); */    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &st );   fmf_freeDestroy( &disG );   return( ret );}

int32 d_piezo_coupling( FMField *out, FMField *strain, FMField *charge_grad,			FMField *mtxG, VolumeGeometry *vg ){  int32 ii, dim, nQP, ret = RET_OK;  FMField *ge = 0, *gptge = 0;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  fmf_createAlloc( &ge, 1, nQP, dim, 1 );  fmf_createAlloc( &gptge, 1, nQP, 1, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( mtxG, ii );    FMF_SetCell( vg->det, ii );    FMF_SetCell( strain, ii );    FMF_SetCell( charge_grad, ii );    fmf_mulAB_nn( ge, mtxG, strain );    fmf_mulATB_nn( gptge, charge_grad, ge );    fmf_sumLevelsMulF( out, gptge, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &ge );  fmf_freeDestroy( &gptge );  return( ret );}

/*!  @par Revision history:  - c: 23.04.2007, r: 23.01.2008*/int32 dw_diffusion_r( FMField *out, FMField *mtxD, Mapping *vg ){  int32 ii, nQP, nEP, ret = RET_OK;  FMField *gtd = 0;  nQP = vg->bfGM->nLev;  nEP = vg->bfGM->nCol;  fmf_createAlloc( &gtd, 1, nQP, nEP, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( vg->bfGM, ii );    FMF_SetCell( vg->det, ii );    if (mtxD->nCell > 1) {      FMF_SetCell( mtxD, ii );    }    fmf_mulATB_nn( gtd, vg->bfGM, mtxD );    fmf_sumLevelsMulF( out, gtd, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &gtd );  return( ret );}

/*!  @par Revision history:  - c: 23.04.2007, r: 23.01.2008*/int32 dw_permeability_r( FMField *out, FMField *mtxD, VolumeGeometry *vg,			 int32 *conn, int32 nEl, int32 nEP,			 int32 *elList, int32 elList_nRow ){  int32 ii, iel, nQP, ret = RET_OK;  FMField *gtd = 0;  nQP = vg->bfGM->nLev;  fmf_createAlloc( &gtd, 1, nQP, nEP, 1 );  for (ii = 0; ii < elList_nRow; ii++) {    iel = elList[ii];    FMF_SetCell( out, ii );    FMF_SetCell( vg->bfGM, iel );    FMF_SetCell( vg->det, iel );    if (mtxD->nCell > 1) {      FMF_SetCell( mtxD, ii );    }    fmf_mulATB_nn( gtd, vg->bfGM, mtxD );    fmf_sumLevelsMulF( out, gtd, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &gtd );   return( ret );}

// `incident` must be preallocated - use mesh_count_incident().int32 mesh_get_incident(Mesh *mesh,                        MeshConnectivity *incident, int32 dim,                        Indices *entities, int32 dent){  int32 ret = RET_OK;  uint32 ii;  uint32 D = mesh->topology->max_dim;  MeshEntityIterator it0[1], it1[1];  MeshConnectivity *conn = mesh->topology->conn[IJ(D, dent, dim)];  if (!conn->num) {    errput("connectivity %d -> %d is not avaliable!/n", dent, dim);    ERR_CheckGo(ret);  }  ii = 0;  incident->offsets[0] = 0;  for (mei_init_sub(it0, mesh, entities, dent); mei_go(it0); mei_next(it0)) {    for (mei_init_conn(it1, it0->entity, dim); mei_go(it1); mei_next(it1)) {      incident->indices[ii++] = it1->entity->ii;    }    incident->offsets[it0->it + 1] = incident->offsets[it0->it] + it1->it_end;  } end_label:  return(ret);}

/*!  @par Revision history:  - c: 12.03.2007, r: 23.01.2008*/int32 d_diffusion( FMField *out, FMField *gradP1, FMField *gradP2,		   FMField *mtxD, Mapping *vg ){  int32 ii, dim, nQP, ret = RET_OK;  FMField *dgp2 = 0, *gp1tdgp2 = 0;  nQP = vg->bfGM->nLev;  dim = vg->bfGM->nRow;  fmf_createAlloc( &dgp2, 1, nQP, dim, 1 );  fmf_createAlloc( &gp1tdgp2, 1, nQP, 1, 1 );  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( vg->det, ii );    FMF_SetCell( gradP1, ii );    FMF_SetCell( gradP2, ii );    if (mtxD->nCell > 1) {      FMF_SetCell( mtxD, ii );    }    fmf_mulAB_nn( dgp2, mtxD, gradP2 );    fmf_mulATB_nn( gp1tdgp2, gradP1, dgp2 );    fmf_sumLevelsMulF( out, gp1tdgp2, vg->det->val );    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &dgp2 );  fmf_freeDestroy( &gp1tdgp2 );  return( ret );}

/*!  For scalar input p: /int_{/Gamma} p dS  For vector input v: /int_{/Gamma} v n dS  @par Revision history:  - 01.11.2007, c*/int32 sg_integrateChunk( SurfaceGeometry *obj, FMField *out, FMField *in,			 int32 *elList, int32 elList_nRow, int32 mode ){  int32 dim, nQP, iel, ii, ret = RET_OK;  FMField *vn = 0;  dim = obj->normal->nRow;  nQP = obj->normal->nLev;     if (in->nRow == 1 || mode == 1) {    for (ii = 0; ii < elList_nRow; ii++) {      iel = elList[ii];      FMF_SetCell( obj->det, iel );      FMF_SetCell( in, ii );      FMF_SetCell( out, ii );      fmf_sumLevelsMulF( out, in, obj->det->val );      ERR_CheckGo( ret );    }  } else if (in->nRow == dim) {    fmf_createAlloc( &vn, 1, nQP, 1, 1 );    for (ii = 0; ii < elList_nRow; ii++) {      iel = elList[ii];      FMF_SetCell( obj->normal, iel );      FMF_SetCell( obj->det, iel );      FMF_SetCell( in, ii );      FMF_SetCell( out, ii );      fmf_mulATB_nn( vn, in, obj->normal );/*       fmf_mulC( vn, -1.0 ); */      fmf_sumLevelsMulF( out, vn, obj->det->val );/*       fmf_print( vn, stdout, 0 ); *//*       fmf_print( out, stdout, 0 ); *//*       sys_pause(); */      ERR_CheckGo( ret );    }  } else {    errput( ErrHead "ERR_Switch/n" );  } end_label:  fmf_freeDestroy( &vn );  return( RET_OK );}

/*!  @par Revision history:  - 24.10.2007, c*/int32 d_sd_st_grad_div( FMField *out, FMField *divU, FMField *gradU,			FMField *divW, FMField *gradW, FMField *divMV,			FMField *gradMV, FMField *coef,			VolumeGeometry *vg_u, int32 mode ){  int32 ii, nQP, ret = RET_OK;  FMField *scalar1 = 0, *scalar2 = 0;  nQP = vg_u->bfGM->nLev;  fmf_createAlloc( &scalar1, 1, nQP, 1, 1 );  if (mode == 1) {    fmf_createAlloc( &scalar2, 1, nQP, 1, 1 );  }  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( coef, ii );    FMF_SetCell( divU, ii );    FMF_SetCell( divW, ii );    FMF_SetCell( vg_u->det, ii );    if (mode == 0) {      // (div u, div w)      fmf_mulATB_nn( scalar1, divU, divW );      fmf_mul( scalar1, coef->val );      fmf_sumLevelsMulF( out, scalar1, vg_u->det->val );    } else if (mode == 1) {      FMF_SetCell( divMV, ii );      FMF_SetCell( gradU, ii );      FMF_SetCell( gradW, ii );      FMF_SetCell( gradMV, ii );      // div u div w div nu ...      fmf_mulATB_nn( scalar1, divU, divW );      fmf_mulATB_nn( scalar2, divMV, scalar1 );      // ... - (gu : gnu^T) div w ...      sub_mul_gradddgrad_scalar( scalar2, gradMV, gradU, divW );      // ... - (gw : gnu^T) div u      sub_mul_gradddgrad_scalar( scalar2, gradMV, gradW, divU );      fmf_mul( scalar2, coef->val );      fmf_sumLevelsMulF( out, scalar2, vg_u->det->val );    }    ERR_CheckGo( ret );  } end_label:  fmf_freeDestroy( &scalar1 );  if (mode == 1) {    fmf_freeDestroy( &scalar2 );  }  return( ret );}

int32 d_diffusion_sa( FMField *out,                      FMField *grad_q, FMField *grad_p,                      FMField *grad_w, FMField *div_w,                      FMField *mtxD, Mapping *vg ){    int32 ii, dim, nQP, ret = RET_OK;    FMField *aux2 = 0, *aux3 = 0, *aux4 = 0, *out0 = 0;    nQP = vg->bfGM->nLev;    dim = vg->bfGM->nRow;    FMF_SetFirst( out );    fmf_createAlloc( &aux2, 1, nQP, dim, 1 );    fmf_createAlloc( &aux3, 1, nQP, 1, 1 );    fmf_createAlloc( &aux4, 1, nQP, dim, 1 );    fmf_createAlloc( &out0, 1, nQP, 1, 1 );    for (ii = 0; ii < out->nCell; ii++) {        FMF_SetCell( vg->bfGM, ii );        FMF_SetCell( vg->det, ii );        FMF_SetCell( mtxD, ii );        FMF_SetCell( grad_q, ii );        FMF_SetCell( grad_p, ii );        FMF_SetCell( grad_w, ii );        FMF_SetCell( div_w, ii );        /* div w K_ij grad_j q grad_i p */        fmf_mulAB_nn( aux2, mtxD, grad_p );        fmf_mulATB_nn( aux3, grad_q, aux2 );        fmf_mulAB_nn( out0, div_w, aux3 );        /* grad_k q K_ij grad_j w_k grad_i p */        fmf_mulATB_nn( aux4, grad_w, aux2 );        fmf_mulATB_nn( aux3, grad_q, aux4 );        fmf_subAB_nn( out0, out0, aux3 );        /* grad_k q K_ij grad_j w_k grad_i p */        fmf_mulAB_nn( aux2, grad_w, grad_p );        fmf_mulAB_nn( aux4, mtxD, aux2 );        fmf_mulATB_nn( aux3, grad_q, aux4 );        fmf_subAB_nn( out0, out0, aux3 );        fmf_sumLevelsMulF( out, out0, vg->det->val );        FMF_SetCellNext( out );        ERR_CheckGo( ret );    }end_label:    fmf_freeDestroy( &out0 );    fmf_freeDestroy( &aux2 );    fmf_freeDestroy( &aux3 );    fmf_freeDestroy( &aux4 );    return( ret );}

int32 mesh_setup_connectivity(Mesh *mesh, int32 d1, int32 d2){  int32 ret = RET_OK;  int32 d3 = 0;  MeshTopology *topology = mesh->topology;  uint32 D = topology->max_dim;  debprintf("request connectivity %d -> %d/n", d1, d2);  if (topology->num[d1] == 0) {    mesh_build(mesh, d1);    ERR_CheckGo(ret);  }  if (topology->num[d2] == 0) {    mesh_build(mesh, d2);    ERR_CheckGo(ret);  }  if (topology->conn[IJ(D, d1, d2)]->num) {    return(ret);  }  if (d1 < d2) {    mesh_setup_connectivity(mesh, d2, d1);    mesh_transpose(mesh, d1, d2);  } else {    if ((d1 == 0) && (d2 == 0)) {      d3 = D;    } else {      d3 = 0;    }    if ((d1 > 0) && (d2 == 0)) {      errput("connectivity %d -> %d should already exist!/n", d1, d2);      ERR_CheckGo(ret);    }    mesh_setup_connectivity(mesh, d1, d3);    mesh_setup_connectivity(mesh, d3, d2);    mesh_intersect(mesh, d1, d2, d3);  }  ERR_CheckGo(ret); end_label:  return(ret);}

/*!  @par Revision history:  - 03.08.2006, c  - 01.12.2006*/int32 dw_biot_grad( FMField *out, float64 coef, FMField *pressure_qp,		    FMField *mtxD, Mapping *svg, Mapping *vvg,		    int32 isDiff ){  int32 ii, nEPU, nEP, dim, nQP, ret = RET_OK;  FMField *dfp = 0, *gtdfp = 0, *gtd = 0, *gtdf = 0;  nQP = vvg->bfGM->nLev;  dim = vvg->bfGM->nRow;  nEPU = vvg->bfGM->nCol;  nEP = svg->bf->nCol;/*   fmf_print( mtxD, stdout, 0 ); */  if (isDiff == 1) {    fmf_createAlloc( &gtd, 1, nQP, dim * nEPU, 1 );    fmf_createAlloc( &gtdf, 1, nQP, dim * nEPU, nEP );  } else {    int32 sym = (dim + 1) * dim / 2;    fmf_createAlloc( &dfp, 1, nQP, sym, 1 );    fmf_createAlloc( &gtdfp, 1, nQP, dim * nEPU, 1 );  }  for (ii = 0; ii < out->nCell; ii++) {    FMF_SetCell( out, ii );    FMF_SetCell( mtxD, ii );    FMF_SetCell( vvg->bfGM, ii );    FMF_SetCell( vvg->det, ii );    if (isDiff == 1) {      FMF_SetCellX1( svg->bf, ii );      form_sdcc_actOpGT_M3( gtd, vvg->bfGM, mtxD );      fmf_mulAB_nn( gtdf, gtd, svg->bf );      fmf_sumLevelsMulF( out, gtdf, vvg->det->val );    } else {      FMF_SetCell( pressure_qp, ii );      fmf_mulAB_nn( dfp, mtxD, pressure_qp );      form_sdcc_actOpGT_VS3( gtdfp, vvg->bfGM, dfp );      fmf_sumLevelsMulF( out, gtdfp, vvg->det->val );    }    ERR_CheckGo( ret );  }  // E.g. 1/dt.  fmfc_mulC( out, coef ); end_label:  if (isDiff) {    fmf_freeDestroy( &gtd );    fmf_freeDestroy( &gtdf );  } else {    fmf_freeDestroy( &dfp );    fmf_freeDestroy( &gtdfp );  }  return( ret );}

// `local_ids` must be preallocated to same size as `entities`.int32 mesh_get_local_ids(Mesh *mesh, Indices *local_ids,                         Indices *entities, int32 dent,                         MeshConnectivity *incident, int32 dim){  int32 ret = RET_OK;  uint32 ii, iind, ic, found;  uint32 D = mesh->topology->max_dim;  MeshEntity entity[1];  MeshEntityIterator it1[1];  MeshConnectivity *conn = mesh->topology->conn[IJ(D, dim, dent)];  if (!conn->num) {    errput("connectivity %d -> %d is not avaliable!/n", dim, dent);    ERR_CheckGo(ret);  }  entity->mesh = mesh;  entity->dim = dim;  ii = 0;  for (iind = 0; iind < incident->num; iind++) {    for (ic = incident->offsets[iind]; ic < incident->offsets[iind+1]; ic++) {      entity->ii = incident->indices[ic];      // printf("%d: ? %d in %d/n", iind, entities->indices[iind], entity->ii);      found = 0;      for (mei_init_conn(it1, entity, dent); mei_go(it1); mei_next(it1)) {        if (entities->indices[iind] == it1->entity->ii) {          local_ids->indices[ii++] = it1->it;          // printf("%d -> %d/n", ii, it1->it);          found = 1;          break; // Degenerate cases - 1. occurrence is returned.        }      }      if (!found) {        errput("entity (%d, %d) not found in entity (%d, %d)!/n",               entities->indices[iind], dent, entity->ii, dim);        ERR_CheckGo(ret);      }    }  } end_label:  return(ret);}

int32 conn_resize(MeshConnectivity *conn, uint32 num, uint32 n_incident){  int32 ret = RET_OK;  conn->num = num;  conn->offsets = realloc_mem(conn->offsets, uint32, num + 1);  ERR_CheckGo(ret);  conn->n_incident = n_incident;  conn->indices = realloc_mem(conn->indices, uint32, n_incident);  ERR_CheckGo(ret); end_label:  if (ERR_Chk) {    errput("conn_resize() failed!");  }  return(ret);}