C++ VecAssemblyBegin函数代码示例

/*---------------------------------------------------------
        Function to fill the rhs vector with
        By + g values ****
---------------------------------------------------------*/
PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
{
  PetscInt       i,j,js,je,jj;
  PetscScalar    val,g[num_z],btri[num_z][3],add_term;
  PetscErrorCode ierr;

  for (i=0; i < nz-2; i++) {
    for (j=0; j <= 2; j++) btri[i][j]=0.0;
    g[i] = 0.0;
  }

  /*  call femBg to set the tri-diagonal b matrix and vector g  */
  femBg(btri,g,nz,z,t);

  /*  setting the entries of the right hand side vector  */
  for (i=0; i < nz-2; i++) {
    val = 0.0;
    js  = 0;
    if (i == 0) js = 1;
    je = 2;
    if (i == nz-2) je = 1;

    for (jj=js; jj <= je; jj++) {
      j    = i+jj-1;
      val += (btri[i][jj])*(y[j]);
    }
    add_term = val + g[i];
    ierr = VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(obj->ksp_rhs);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(obj->ksp_rhs);CHKERRQ(ierr);

  /*  return to main driver function  */
  return 0;
}

PetscErrorCode BSSPotR1Vec(BSS self, PF pot, Vec V) {
  int order = self->order;
  int nb = self->num_basis;
  int ne = self->num_ele;
  int nq = order*ne;

  PetscErrorCode ierr;
  ierr = BSSCheck(self); CHKERRQ(ierr);

  PetscScalar *vs; PetscMalloc1(nq, &vs);  
  ierr = PFApply(pot, nq, self->Rrs, vs); CHKERRQ(ierr);
  for(int i = 0; i < nb; i++) {
    int k0, k1; Non0QuadIndex(i, i, order, nq, &k0, &k1);
    PetscScalar v = 0.0;
    for(int k = k0; k < k1; k++) {
      //v += self->vals[k+i*nq] * self->ws[k] * vs[k];
      v += self->vals[k+i*nq] * self->ws[k] * self->qrs[k] * vs[k];
    }
    ierr = VecSetValue(V, i, v, INSERT_VALUES); CHKERRQ(ierr);
  }

  PetscFree(vs);
  VecAssemblyBegin(V);
  VecAssemblyEnd(V);

  return 0;
}

Ensure(FFT, transforms_constant_into_delta_function)
{
    Vec v;

    DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,2,1,NULL,&da);

    DMCreateGlobalVector(da, &v);
    VecZeroEntries(v);
    VecSetValue(v, 0, 0.7, INSERT_VALUES);
    VecAssemblyBegin(v);
    VecAssemblyEnd(v);

    scFftCreate(da, &fft);

    Vec x, y, z;
    scFftCreateVecsFFTW(fft, &x, &y, &z);

    scFftTransform(fft, v, 0, y);
    const PetscScalar *arr;
    VecGetArrayRead(y, &arr);
    //assert_that_double(fabs(arr[0] - 0.7), is_less_than_double(1.0e-6));
    assert_that(fabs(arr[1] - 0.0) < 1.0e-6, is_true);
    assert_that(fabs(arr[1] - 0.7) < 1.0e-6, is_true);
    assert_that(fabs(arr[2] - 0.0) < 1.0e-6, is_true);
    assert_that(fabs(arr[3] - 0.7) < 1.0e-6, is_true);
    assert_that(fabs(arr[4] - 0.0) < 1.0e-6, is_true);
    VecRestoreArrayRead(y, &arr);

    VecDestroy(&x);
    VecDestroy(&y);
    VecDestroy(&z);
    VecDestroy(&v);
    scFftDestroy(&fft);
    DMDestroy(&da);
}

void operator| (PUP::er& p, Vec& v) {
  PetscInt sz;
  if (!p.isUnpacking()) {
    VecGetSize(v, &sz);
  }
  p | sz;
  if(p.isUnpacking()) {    
	  VecCreateSeq(PETSC_COMM_WORLD, sz,  &v);
	  VecSetFromOptions(v);
	  VecGetSize(v, &sz);

	  for (int i = 0; i < sz; i++) {
		  PetscScalar d;
		  p | d;
		  VecSetValue(v, i, d, INSERT_VALUES);
	  }
	  VecAssemblyBegin(v);
	  VecAssemblyEnd(v);
  }
  else {
	  for (int i = 0; i < sz; i++) {
		  PetscScalar d;
		  VecGetValues(v, 1, &i, &d);
		  p | d;
	  }
  }
}

double Assembly_Mat_Vec(InterpolationDataStruct* pIData, int field, Mat M, Vec Fx, Vec Fy) {
	double coord1[3],coord2[3],coord3[3], value1, value2, value3, Fx_values[3],Fy_values[3];

	//Assembly
	FIter fit = M_faceIter(pIData->m2);
	while (pEntity entity = FIter_next(fit)){

		//get the coordinates and the id's from the element's vertices
		pEntity vertices[3] = {entity->get(0,0), entity->get(0,1), entity->get(0,2)};
		int IDs[3] = {EN_id(vertices[0]), EN_id(vertices[1]), EN_id(vertices[2])};
		int position[3] = {IDs[0]-1, IDs[1]-1, IDs[2]-1};
		V_coord(vertices[0],coord1);
		V_coord(vertices[1],coord2);
		V_coord(vertices[2],coord3);

		//calculate area
		double At = F_area(coord1,coord2,coord3);
		double At_by_6 = (double)(At/6.);
		double At_by_12 = (double)(At/12.);
		double M_values[9] = {At_by_6, At_by_12, At_by_12, At_by_12, At_by_6, At_by_12, At_by_12, At_by_12, At_by_6};

		//get node values
		pIData->pGetDblFunctions[field](IDs[0]-1,value1);
		pIData->pGetDblFunctions[field](IDs[1]-1,value2);
		pIData->pGetDblFunctions[field](IDs[2]-1,value3);
		double CV[3] = {value1/6., value2/6., value3/6.};

		//calculate the contributions
		for(int i=0;i<3;i++){
			Fx_values[i] = (coord2[1]-coord3[1])*CV[0] + (coord3[1]-coord1[1])*CV[1] + (coord1[1]-coord2[1])*CV[2];
			Fy_values[i] = (coord3[0]-coord2[0])*CV[0] + (coord1[0]-coord3[0])*CV[1] + (coord2[0]-coord1[0])*CV[2];
		}
		VecSetValues(Fx,3,position,Fx_values, ADD_VALUES);
		VecSetValues(Fy,3,position,Fy_values, ADD_VALUES);
		MatSetValues(M,3,position,3,position,M_values, ADD_VALUES);
	}
	FIter_delete(fit);

	PetscErrorCode ierr;
	ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	ierr = VecAssemblyBegin(Fx);CHKERRQ(ierr);
	ierr = VecAssemblyEnd(Fx);CHKERRQ(ierr);
	ierr = VecAssemblyBegin(Fy);CHKERRQ(ierr);
	ierr = VecAssemblyEnd(Fy);CHKERRQ(ierr);
	return 0;
}

/*@
    KSPComputeExplicitOperator - Computes the explicit preconditioned operator.

    Collective on KSP

    Input Parameter:
.   ksp - the Krylov subspace context

    Output Parameter:
.   mat - the explict preconditioned operator

    Notes:
    This computation is done by applying the operators to columns of the
    identity matrix.

    Currently, this routine uses a dense matrix format when 1 processor
    is used and a sparse format otherwise.  This routine is costly in general,
    and is recommended for use only with relatively small systems.

    Level: advanced

.keywords: KSP, compute, explicit, operator

.seealso: KSPComputeEigenvaluesExplicitly(), PCComputeExplicitOperator()
@*/
PetscErrorCode  KSPComputeExplicitOperator(KSP ksp,Mat *mat)
{
  Vec            in,out;
  PetscErrorCode ierr;
  PetscMPIInt    size;
  PetscInt       i,M,m,*rows,start,end;
  Mat            A;
  MPI_Comm       comm;
  PetscScalar    *array,one = 1.0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
  PetscValidPointer(mat,2);
  ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = VecDuplicate(ksp->vec_sol,&in);CHKERRQ(ierr);
  ierr = VecDuplicate(ksp->vec_sol,&out);CHKERRQ(ierr);
  ierr = VecGetSize(in,&M);CHKERRQ(ierr);
  ierr = VecGetLocalSize(in,&m);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(in,&start,&end);CHKERRQ(ierr);
  ierr = PetscMalloc1(m,&rows);CHKERRQ(ierr);
  for (i=0; i<m; i++) rows[i] = start + i;

  ierr = MatCreate(comm,mat);CHKERRQ(ierr);
  ierr = MatSetSizes(*mat,m,m,M,M);CHKERRQ(ierr);
  if (size == 1) {
    ierr = MatSetType(*mat,MATSEQDENSE);CHKERRQ(ierr);
    ierr = MatSeqDenseSetPreallocation(*mat,NULL);CHKERRQ(ierr);
  } else {
    ierr = MatSetType(*mat,MATMPIAIJ);CHKERRQ(ierr);
    ierr = MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);CHKERRQ(ierr);
  }
  ierr = MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
  ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr);

  for (i=0; i<M; i++) {

    ierr = VecSet(in,0.0);CHKERRQ(ierr);
    ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(in);CHKERRQ(ierr);

    ierr = KSP_MatMult(ksp,A,in,out);CHKERRQ(ierr);
    ierr = KSP_PCApply(ksp,out,in);CHKERRQ(ierr);

    ierr = VecGetArray(in,&array);CHKERRQ(ierr);
    ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecRestoreArray(in,&array);CHKERRQ(ierr);

  }
  ierr = PetscFree(rows);CHKERRQ(ierr);
  ierr = VecDestroy(&in);CHKERRQ(ierr);
  ierr = VecDestroy(&out);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

double linearSystemPETSc<scalar>::normInfRightHandSide() const
{
  PetscReal nor;
  VecAssemblyBegin(_b);
  VecAssemblyEnd(_b);
  _try(VecNorm(_b, NORM_INFINITY, &nor));
  return nor;
}

void linearSystemPETSc<scalar>::zeroRightHandSide()
{
  if (_isAllocated) {
    _try(VecAssemblyBegin(_b));
    _try(VecAssemblyEnd(_b));
    _try(VecZeroEntries(_b));
  }
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscInt       i,blocks[2],nlocal;
  PetscMPIInt    size,rank;
  PetscScalar    value;
  Vec            x,y;
  IS             is1,is2;
  VecScatter     ctx = 0;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (size != 2) SETERRQ(PETSC_COMM_SELF,1,"Must run with 2 processors");

  /* create two vectors */
  if (!rank) nlocal = 8;
  else nlocal = 4;
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,nlocal,12);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,8,&y);CHKERRQ(ierr);

  /* create two index sets */
  if (!rank) {
    blocks[0] = 0; blocks[1] = 2;
  } else {
    blocks[0] = 1; blocks[1] = 2;
  }
  ierr = ISCreateBlock(PETSC_COMM_SELF,4,2,blocks,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,8,0,1,&is2);CHKERRQ(ierr);

  for (i=0; i<12; i++) {
    value = i;
    ierr = VecSetValues(x,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x);CHKERRQ(ierr);

  ierr = VecScatterCreate(x,is1,y,is2,&ctx);CHKERRQ(ierr);
  ierr = VecScatterBegin(ctx,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(ctx,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr);

  ierr = PetscSleep(2.0*rank);CHKERRQ(ierr);
  ierr = VecView(y,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);

  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&y);CHKERRQ(ierr);
  ierr = ISDestroy(&is1);CHKERRQ(ierr);
  ierr = ISDestroy(&is2);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return 0;
}

/*
#undef __FUNCT__
#define __FUNCT__ "phi2eval"
PetscErrorCode phi2eval(Vec *X, Vec *Phi, PetscInt n) {
  PetscInt i,j,k,lo,hi;
  PetscErrorCode ierr;
  PetscReal sqrt2;

  PetscFunctionBegin;
  sqrt = PetscSqrtReal(2.0);
  ierr = VecGetOwnershipRange(X,&lo,&hi);CHKERRQ(ierr);
PetscErrorCode phi2eval(PetscReal *x, PetscInt n, PetscReal *phi) {
// Phi = .5*[x(1)^2  sqrt(2)*x(1)*x(2) ... sqrt(2)*x(1)*x(n) ... x(2)^2 sqrt(2)*x(2)*x(3) .. x(n)^2]
    PetscInt i,j,k;
    PetscReal sqrt2 = PetscSqrtReal(2.0);
    PetscFunctionBegin;
    j=0;

    for (i=0;i<n;i++) {
        phi[j] = 0.5 * x[i]*x[i];
        j++;
        for (k=i+1;k<n;k++) {
            phi[j]  = x[i]*x[k]/sqrt2;
            j++;
        }

    }

    PetscFunctionReturn(0);
    }*/
PetscErrorCode PoundersGramSchmidtReset(TAO_POUNDERS *mfqP, Vec *Q, PetscInt n)
{
  PetscInt i;
  for (i=0;i<n;i++) {
    ierr = VecSet(Q[i],0.0);CHKERRQ(ierr);
    ierr = VecSetValue(Q[i],i,1.0,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(Q[i]);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(Q[i]);CHKERRQ(ierr);
  }
}

   -fout <file> : output file name\n\n";

#include <petscmat.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Mat            A;
  Vec            b;
  char           fileout[PETSC_MAX_PATH_LEN];
  PetscInt       i,j,m = 6,n = 6,N = 36,Ii,J;
  PetscErrorCode ierr;
  PetscScalar    val,v;
  PetscViewer    view;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      v = -1.0;  Ii = j + n*i;
      if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
      if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
      if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
      if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
      v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
    }
  }
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecSetSizes(b,PETSC_DECIDE,N);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
  for (i=0; i<N; i++) {
    val  = i + 1;
    ierr = VecSetValues(b,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

  ierr = PetscOptionsGetString(NULL,NULL,"-fout",fileout,PETSC_MAX_PATH_LEN,NULL);CHKERRQ(ierr);
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileout,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
  ierr = MatView(A,view);CHKERRQ(ierr);
  ierr = VecView(b,view);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);

  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}

/*@
    MatComputeExplicitOperator - Computes the explicit matrix

    Collective on Mat

    Input Parameter:
.   inmat - the matrix

    Output Parameter:
.   mat - the explict preconditioned operator

    Notes:
    This computation is done by applying the operators to columns of the
    identity matrix.

    Currently, this routine uses a dense matrix format when 1 processor
    is used and a sparse format otherwise.  This routine is costly in general,
    and is recommended for use only with relatively small systems.

    Level: advanced

.keywords: Mat, compute, explicit, operator
@*/
PetscErrorCode  MatComputeExplicitOperator(Mat inmat,Mat *mat)
{
  Vec            in,out;
  PetscErrorCode ierr;
  PetscInt       i,m,n,M,N,*rows,start,end;
  MPI_Comm       comm;
  PetscScalar    *array,zero = 0.0,one = 1.0;
  PetscMPIInt    size;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(inmat,MAT_CLASSID,1);
  PetscValidPointer(mat,2);

  ierr = PetscObjectGetComm((PetscObject)inmat,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = MatGetLocalSize(inmat,&m,&n);CHKERRQ(ierr);
  ierr = MatGetSize(inmat,&M,&N);CHKERRQ(ierr);
  ierr = MatGetVecs(inmat,&in,&out);CHKERRQ(ierr);
  ierr = VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(out,&start,&end);CHKERRQ(ierr);
  ierr = PetscMalloc(m*sizeof(PetscInt),&rows);CHKERRQ(ierr);
  for (i=0; i<m; i++) rows[i] = start + i;

  ierr = MatCreate(comm,mat);CHKERRQ(ierr);
  ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr);
  if (size == 1) {
    ierr = MatSetType(*mat,MATSEQDENSE);CHKERRQ(ierr);
    ierr = MatSeqDenseSetPreallocation(*mat,NULL);CHKERRQ(ierr);
  } else {
    ierr = MatSetType(*mat,MATMPIAIJ);CHKERRQ(ierr);
    ierr = MatMPIAIJSetPreallocation(*mat,n,NULL,N-n,NULL);CHKERRQ(ierr);
  }

  for (i=0; i<N; i++) {

    ierr = VecSet(in,zero);CHKERRQ(ierr);
    ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(in);CHKERRQ(ierr);

    ierr = MatMult(inmat,in,out);CHKERRQ(ierr);

    ierr = VecGetArray(out,&array);CHKERRQ(ierr);
    ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecRestoreArray(out,&array);CHKERRQ(ierr);

  }
  ierr = PetscFree(rows);CHKERRQ(ierr);
  ierr = VecDestroy(&out);CHKERRQ(ierr);
  ierr = VecDestroy(&in);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

void FEMAssemble2DElasticity(MPI_Comm comm, Mesh *mesh, Mat &A, Vec &b,
                             PetscReal E, PetscReal mi, PetscReal(*dens)(Element *),
                             void(*f)(Element*, PetscReal, PetscReal*)) {



    PetscInt rank;
    MPI_Comm_rank(comm, &rank);

    PetscInt size = mesh->vetrices.size();


    MatCreateSeqAIJ(PETSC_COMM_SELF, size * 2, size * 2, 14, PETSC_NULL, &A);
    VecCreateMPI(comm, size * 2, PETSC_DECIDE, &b);



    for (std::map<PetscInt, Element*>::iterator e = mesh->elements.begin(); e
            != mesh->elements.end(); e++) {



        Point* vetrices[3];
        PetscInt vIndex[3];
        for (int i = 0; i < 3; i++) {
            vIndex[i] = e->second->vetrices[i];
            vetrices[i] = mesh->vetrices[vIndex[i]];
        }
        PetscReal lStiff[36];
        PetscReal bLoc[6];

        //Local stiffness matrix and force vector assembly
        PetscReal fs[2];
        f(e->second, dens(e->second), fs);
        elastLoc(vetrices, E, mi, fs, lStiff, bLoc);

        PetscInt idx[6], idxLoc[6];
        for (int i = 0; i < 3; i++) {
            idx[i * 2] = vIndex[i] * 2;
            idx[i * 2 + 1] = vIndex[i] * 2 + 1;

            idxLoc[i * 2] = vIndex[i] * 2 - mesh->startIndexes[rank] * 2;
            idxLoc[i * 2 + 1] = vIndex[i] * 2 + 1 - mesh->startIndexes[rank] * 2;
        }

        MatSetValues(A, 6, idxLoc, 6, idxLoc, lStiff, ADD_VALUES);
        VecSetValues(b, 6, idx, bLoc, ADD_VALUES);
    }

    VecAssemblyBegin(b);
    VecAssemblyEnd(b);

    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
}

int main(int argc,char **args)
{
  Mat            mat;
  Vec            b,u;
  PC             pc;
  PetscErrorCode ierr;
  PetscInt       n = 5,i,col[3];
  PetscScalar    value[3];

  PetscInitialize(&argc,&args,(char *)0,help);

  /* Create vectors */
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&b);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&u);CHKERRQ(ierr);

  /* Create and assemble matrix */
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,n,n,PETSC_NULL,&mat);CHKERRQ(ierr);
  value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
  for (i=1; i<n-1; i++) {
    col[0] = i-1; col[1] = i; col[2] = i+1;
    ierr = MatSetValues(mat,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  }
  i = n - 1; col[0] = n - 2; col[1] = n - 1;
  ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
  ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* Create PC context and set up data structures */
  ierr = PCCreate(PETSC_COMM_WORLD,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCSOR);CHKERRQ(ierr);
  ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
  ierr = PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = PCSetUp(pc);CHKERRQ(ierr);

  value[0] = 1.0;
  for (i=0; i<n; i++) {
    ierr = VecSet(u,0.0);CHKERRQ(ierr);
    ierr = VecSetValues(u,1,&i,value,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
    ierr = PCApply(pc,u,b);CHKERRQ(ierr);
    ierr = VecView(b,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  }

  /* Free data structures */
  ierr = MatDestroy(&mat);CHKERRQ(ierr);
  ierr = PCDestroy(&pc);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr); 
  ierr = PetscFinalize();
  return 0;
}

PetscErrorCode FEMAssembleTotal2DLaplace(MPI_Comm comm, Mesh *mesh, Mat &A,
        Vec &b, PetscReal(*f)(Point), PetscReal(*K)(Point)) {
    PetscErrorCode ierr;
    PetscInt rank;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

    PetscInt size = mesh->vetrices.size();

    ierr
        = MatCreateMPIAIJ(comm, size, size, PETSC_DECIDE, PETSC_DECIDE, 7, PETSC_NULL, 0, PETSC_NULL, &A);
    CHKERRQ(ierr);
    ierr = VecCreateMPI(comm, size, PETSC_DECIDE, &b);
    CHKERRQ(ierr);

    for (std::map<PetscInt, Element*>::iterator e = mesh->elements.begin(); e
            != mesh->elements.end(); e++) {
        PetscScalar bl[3];
        PetscScalar Al[9];
        PetscReal R[4];

        PetscInt elSize = e->second->numVetrices;
        Point *vetrices = new Point[elSize];
        PetscInt ixs[elSize];
        for (int j = 0; j < elSize; j++) {
            ixs[j] = e->second->vetrices[j];
            vetrices[j] = *(mesh->vetrices[ixs[j]]);
        }

        R[0] = vetrices[1].x - vetrices[0].x;
        R[2] = vetrices[1].y - vetrices[0].y;
        R[1] = vetrices[2].x - vetrices[0].x;
        R[3] = vetrices[2].y - vetrices[0].y;

        Point center = getCenterOfSet(vetrices, elSize);
        bLoc(R, bl, f(center));
        ALoc(R, Al, K(center));

        ierr = VecSetValues(b, elSize, ixs, bl, ADD_VALUES);
        CHKERRQ(ierr);
        ierr = MatSetValues(A, elSize, ixs, elSize, ixs, Al, ADD_VALUES);
        CHKERRQ(ierr);
    }
    ierr = VecAssemblyBegin(b);
    CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);
    CHKERRQ(ierr);

    ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    return ierr;

}