facnb.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: facnb.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "init.h"
#include "util/util.h"

#define USEQSORT
#ifdef USEQSORT
#include <stdlib.h>

int *facvert;

int CompareFaces(const void* v1,const void* v2) {
  const int* t1=(const int*)v1;
  const int* t2=(const int*)v2;

  if (facvert[NN* *t1+0]<facvert[NN* *t2+0]) return -1;
  if (facvert[NN* *t1+0]>facvert[NN* *t2+0]) return 1;

  if (facvert[NN* *t1+1]<facvert[NN* *t2+1]) return -1;
  if (facvert[NN* *t1+1]>facvert[NN* *t2+1]) return -1;

  if (facvert[NN* *t1+2]<facvert[NN* *t2+2]) return -1;
  if (facvert[NN* *t1+2]>facvert[NN* *t2+2]) return -1;

  /* we should never reach this point */
  assert(0==1);
  return 0;
}
#endif

int FacNB(GridData *gdata)
{
  int  t_facvert[NF*NV];
  int  t_elevs[NV];

  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (rank) {
    gdata->ln_bnd_fac=0;
    MagparFunctionLogReturn(0);
  }

  /* get connectivity matrix of elements */
  /* recommended size (METIS manual):
     sizeof(ia): (n_ele+1)
     sizeof(ja): 4*n_ele
  */
  int *ia,*ja;
  ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);

  for (int i=1;i<gdata->n_ele+1;i++) {
    if (ia[i]-ia[i-1]>NF) {
      SETERRQ2(PETSC_ERR_MEMC,"Found element with %i neighbors (>NF)!\n",ia[i]-ia[i-1],NF);
    }
  }
  assert(ia[gdata->n_ele]<NF*gdata->n_ele);

  /* check for array overrun */
  assert(ia[gdata->n_ele]<NF*gdata->n_ele);

  int *elefac;
  ierr = PetscMalloc(gdata->n_ele*NF*sizeof(int),&elefac);CHKERRQ(ierr);

  int *facnb;
  ierr = PetscMalloc(NF*2*gdata->n_ele*sizeof(int),&facnb);CHKERRQ(ierr);

  for (int i=0;i<NF*2*gdata->n_ele;i++) {
    facnb[i]=C_UNK;
  }

  /* each element has NF faces,
     the last entry in ia (ia[gdata->n_ele])
     gives us the total number of neighbours,
     each pair of neighbours appears twice,
     thus we have ia[gdata->n_ele]/2 shared faces,
     thus the total number of faces is
     NF*n_ele-ia[gdata->n_ele]/2
  */
  int n_fac;
  n_fac=NF*gdata->n_ele-ia[gdata->n_ele]/2;

#ifndef USEQSORT
  int *facvert;
#endif
  ierr = PetscMalloc(n_fac*NN*sizeof(int),&facvert);CHKERRQ(ierr);

  int t_facmax;
  t_facmax=0;

  ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Searching for surface faces and nodes...\n");
  for (int i=0;i<gdata->n_ele;i++) {
    ierr = ProgressBar(i,gdata->n_ele,10);
    /* get number of neighbours for current element */
    int t_nnb;
    t_nnb=ia[i+1]-ia[i];

    /* get nodes */
    for (int j=0; j<NV; j++) {
      t_elevs[j]=gdata->elevert[NV*i+j];
/*
      if (mergepoints) {
        t_elevs[j]=PetscMin(gdata->elevert[NV*i+j],gdata->vertcoupl[gdata->elevert[NV*i+j]]);
      }
*/
    }

    /* sort nodes */
    ierr = PetscSortInt(NV,(PetscInt*)t_elevs);CHKERRQ(ierr);

    /* get facvert */
    t_facvert[0]  = t_elevs[0];  /* first face */
    t_facvert[1]  = t_elevs[1];
    t_facvert[2]  = t_elevs[2];
    t_facvert[3]  = t_elevs[3];

    t_facvert[4]  = t_elevs[0];  /* second face */
    t_facvert[5]  = t_elevs[1];
    t_facvert[6]  = t_elevs[3];
    t_facvert[7]  = t_elevs[2];

    t_facvert[8]  = t_elevs[0];  /* third face */
    t_facvert[9]  = t_elevs[2];
    t_facvert[10] = t_elevs[3];
    t_facvert[11] = t_elevs[1];

    t_facvert[12] = t_elevs[1];  /* fourth face */
    t_facvert[13] = t_elevs[2];
    t_facvert[14] = t_elevs[3];
    t_facvert[15] = t_elevs[0];

    /* process all NF faces */
    for (int j=0; j<NF; j++) {
      int t_found;
      t_found=0;

      /* search the neighbours */
      for (int k=0; k<t_nnb; k++) {
        /* get the id of our neighbour */
        int t_nbid;
        t_nbid=ja[ia[i]+k];

        /* find shared nodes */
        int x,vshared[NV];
        x=tetnb(gdata->elevert+NV*i,gdata->elevert+NV*t_nbid,vshared);
        /* we assume we get neighbors sharing a face - i.e. 3 nodes - from Metis */
        assert(x==NN);

        if ((t_facvert[NV*j+0]==vshared[0]) &&
            (t_facvert[NV*j+1]==vshared[1]) &&
            (t_facvert[NV*j+2]==vshared[2])
           ) {

           t_found=1;

          /* if the id of the neighbour is smaller than ours
             then its faces have already been processed */
          if (t_nbid<i) {
            /* search for our id among neighbours of t_nbid */
            for (int m=0; m<NF; m++) {
              if (facnb[elefac[t_nbid*NF+m]*2+1]==i) {
                elefac[i*NF+j]=elefac[t_nbid*NF+m];
                break;
              }
            }
          }
          /* else we have to register a new face */
          else {
            facvert[t_facmax*NN+0]=vshared[0];
            facvert[t_facmax*NN+1]=vshared[1];
            facvert[t_facmax*NN+2]=vshared[2];
            facnb[t_facmax*2+0]=i;
            facnb[t_facmax*2+1]=t_nbid;
            elefac[i*NF+j]=t_facmax;
            t_facmax++;
          }
          break;
        }
        if (t_found) break;
      }
      /* if we cannot find a neighbour, it must be on the surface */
      if (!t_found) {
        /* check orientation (normal vector should point outwards),
           this ensures, that surface triangles have the correct
           orientation
        */
        /* if a surface triangle is found the element must have
           less than NF neighbors
        */
        assert(t_nnb<NF);

        PetscReal *x1,*x2,*x3,*x4;
        x1=gdata->vertxyz+(t_facvert[NV*j+0])*ND;
        x2=gdata->vertxyz+(t_facvert[NV*j+1])*ND;
        x3=gdata->vertxyz+(t_facvert[NV*j+2])*ND;
        x4=gdata->vertxyz+(t_facvert[NV*j+3])*ND;

        facvert[t_facmax*NN+0]=t_facvert[NV*j+0];
        /* proper orientation */
        if (tetvol(x1,x2,x3,x4)<0.0) {
          facvert[t_facmax*NN+1]=t_facvert[NV*j+1];
          facvert[t_facmax*NN+2]=t_facvert[NV*j+2];
        }
        /* swap second and third index vertex id */
        else {
          facvert[t_facmax*NN+1]=t_facvert[NV*j+2];
          facvert[t_facmax*NN+2]=t_facvert[NV*j+1];
        }
        facnb[t_facmax*2+0]=i;
        facnb[t_facmax*2+1]=C_BND;
        elefac[i*NF+j]=t_facmax;
        t_facmax++;
      }
    }
  }
  ierr = ProgressBar(1,1,10);

  PetscLogDouble t_t3;
  ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "Search for surface faces and nodes took %g s = %g min\n",
    t_t3-t_t2,
    (t_t3-t_t2)/60.0
  );

#ifdef ADDONS
  PetscTruth flg;
  int mergepoints;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-hdemag_mergepoints",(PetscInt*)&mergepoints,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    mergepoints=1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hdemag_mergepoints not found, using default value: %i\n",
      mergepoints
    );
  }
  if (gdata->vertcoupl==PETSC_NULL) {
    mergepoints=0;
    PetscPrintf(PETSC_COMM_WORLD,"No coupled vertices found, setting -hdemag_mergepoints %i!\n",mergepoints);
  }

  int t_cnt;
  t_cnt=0;
  if (mergepoints) {
    /* update facvert array with coupled nodes */
    for (int i=0;i<n_fac;i++) {
      for (int j=0;j<NN;j++) {
        if (gdata->vertcoupl[facvert[NN*i+j]]>=0) {
          /* use main coupled node (lowest id) */
          facvert[NN*i+j]=PetscMin(facvert[NN*i+j],gdata->vertcoupl[facvert[NN*i+j]]);
        }
      }
#define swap(a,b) {int tswap;tswap=a;a=b;b=tswap;}
      if (facvert[NN*i+0]>facvert[NN*i+1]) {
        swap(facvert[NN*i+0],facvert[NN*i+1]);
        swap(facvert[NN*i+1],facvert[NN*i+2]);
      }
      if (facvert[NN*i+0]>facvert[NN*i+2]) {
        swap(facvert[NN*i+0],facvert[NN*i+1]);
        swap(facvert[NN*i+1],facvert[NN*i+2]);
      }
      if (facvert[NN*i+0]>facvert[NN*i+1]) {
        swap(facvert[NN*i+0],facvert[NN*i+1]);
        swap(facvert[NN*i+1],facvert[NN*i+2]);
      }
    }

    PetscPrintf(PETSC_COMM_WORLD,"Searching for coupled faces/nodes...\n");

    int nsfac=0;
    int *sfac;
    ierr = PetscMalloc(n_fac*sizeof(int),&sfac);CHKERRQ(ierr);

    for (int i=0;i<n_fac;i++) {
      if (facnb[i*2+1]!=C_BND ||
          gdata->vertcoupl[facvert[NN*i+0]]==-1 ||
          gdata->vertcoupl[facvert[NN*i+1]]==-1 ||
          gdata->vertcoupl[facvert[NN*i+2]]==-1) continue;

      sfac[nsfac]=i;
      nsfac++;
    }

#ifdef USEQSORT
    PetscPrintf(PETSC_COMM_WORLD,"Sorting surface triangles...\n");
    qsort(sfac,nsfac,sizeof(int),CompareFaces);
#endif

    PetscPrintf(PETSC_COMM_WORLD,"Checking %i coupled faces/nodes...\n",nsfac);

    for (int i2=0;i2<nsfac-1;i2++) {
      ierr = ProgressBar(i2,nsfac-1,10);

      int i;
      i=sfac[i2];

      /* search all other boundary faces */
      for (int k2=i2+1;k2<nsfac;k2++) {

        int k;
        k=sfac[k2];

#ifdef USEQSORT
        if (facvert[NN*k+0]>facvert[NN*i+0]) break;
#endif
        /* if all three vertices are coupled */
        /* we need to check only three cases, because
           faces are oriented to give positive volume
           and our coupled faces must be oppositely oriented
           and vertices are sorted
        */
        if (
          facvert[NN*i+0]==facvert[NN*k+0] &&
          facvert[NN*i+1]==facvert[NN*k+2] &&
          facvert[NN*i+2]==facvert[NN*k+1]
          ) {
          /* mark faces as interior but coupled */
          facnb[i*2+1]=facnb[k*2+0];
          facnb[k*2+1]=facnb[i*2+0];

          t_cnt+=2;
          break;
        }
      }
    }
    ierr = ProgressBar(1,1,10);

    ierr = PetscFree(sfac);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"Retagged %i coupled surfaces\n",t_cnt);
  }

  t_t2=t_t3;
  ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "Search for coupled faces/nodes took %g s = %g min\n",
    t_t3-t_t2,
    (t_t3-t_t2)/60.0
  );
#endif

  ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertbndg2bnd);CHKERRQ(ierr);

  /* reset boundary indicators */
  for (int i=0;i<gdata->n_vert;i++) {
    gdata->vertbndg2bnd[i]=C_INT;
  }

  /* set boundary indicators and count boundary vertices */
  gdata->n_vert_bnd=0;
  for (int i=0;i<n_fac;i++) {
    if (facnb[i*2+1]==C_BND) {
      for (int j=0;j<NN;j++) {
        int k;
        k=facvert[NN*i+j];
#ifdef ADDONS
        if (mergepoints && gdata->vertcoupl[k]>=0) {
          /* check that k is the main coupling node (lowest id) as set above */
          assert(k<gdata->vertcoupl[k]);
        }
#endif
        /* avoid counting vertices twice */
        if (gdata->vertbndg2bnd[k]==C_INT) {
          gdata->vertbndg2bnd[k]=C_BND;
          gdata->n_vert_bnd++;
        }
      }
    }
  }
  PetscPrintf(PETSC_COMM_WORLD,
    "Boundary matrix: %i vertices, %g matrix elements, %i MB \n",
    gdata->n_vert_bnd,
    1.0*gdata->n_vert_bnd*gdata->n_vert_bnd,
    gdata->n_vert_bnd/1024*gdata->n_vert_bnd/1024*sizeof(PetscReal)+1
  );

  /* count boundary faces */
  gdata->n_bnd_fac=0;
  for (int i=0;i<gdata->n_ele;i++) {
    for (int j=0;j<NF;j++) {
      if (facnb[elefac[i*NF+j]*2+1]==C_BND) {
        gdata->n_bnd_fac++;
      }
    }
  }
  gdata->ln_bnd_fac=gdata->n_bnd_fac;
  PetscPrintf(PETSC_COMM_WORLD,"Number of boundary faces: %i\n",gdata->n_bnd_fac);

  /* store data of boundary faces */
  ierr = PetscMalloc(gdata->n_bnd_fac*NN*sizeof(int),&gdata->bndfacvert);CHKERRQ(ierr);
  int t_cnt2;
  t_cnt2=0;
  for (int i=0;i<n_fac;i++) {
    if (facnb[i*2+1]==C_BND) {
      for (int j=0;j<NN;j++) {
        gdata->bndfacvert[t_cnt2*NN+j]=facvert[i*NN+j];
      }
      t_cnt2++;
    }
  }
  assert(t_cnt2==gdata->n_bnd_fac);

  /* convert boundary indicators to new boundary node numbers */
  int t_cnt3;
  t_cnt3=0;
  for (int i=0; i<gdata->n_vert; i++) {
    if (gdata->vertbndg2bnd[i]==C_BND) {
      gdata->vertbndg2bnd[i]=t_cnt3++;
    }
  }
  assert(t_cnt3==gdata->n_vert_bnd);

  ierr = PetscFree(elefac);CHKERRQ(ierr);
  ierr = PetscFree(facnb);CHKERRQ(ierr);
  ierr = PetscFree(facvert);CHKERRQ(ierr);

  ierr = PetscFree(ia);CHKERRQ(ierr);
  ierr = PetscFree(ja);CHKERRQ(ierr);


  MagparFunctionLogReturn(0);
}