Constrain Projection Matrix

Collaboration diagram for Constrain Projection Matrix:

Classes
struct	ConstrainMatrixCtx
	structure for projection matrices More...

Functions
MoFEMErrorCode	ProjectionMatrixMultOpQ (Mat Q, Vec x, Vec f)
	Multiplication operator for Q = I-CTC(CCT)^-1C. More...
MoFEMErrorCode	ConstrainMatrixMultOpP (Mat P, Vec x, Vec f)
	Multiplication operator for P = CT(CCT)^-1C. More...
MoFEMErrorCode	ConstrainMatrixMultOpR (Mat R, Vec x, Vec f)
	Multiplication operator for R = CT(CCT)^-1. More...
MoFEMErrorCode	ConstrainMatrixMultOpRT (Mat RT, Vec x, Vec f)
	Multiplication operator for RT = (CCT)^-TC. More...
MoFEMErrorCode	ConstrainMatrixMultOpCTC_QTKQ (Mat CTC_QTKQ, Vec x, Vec f)
	Multiplication operator for RT = (CCT)^-TC. More...
MoFEMErrorCode	ConstrainMatrixDestroyOpPorQ (Mat Q)
	Destroy shell matrix Q. More...
MoFEMErrorCode	ConstrainMatrixDestroyOpQTKQ (Mat QTKQ)
	Destroy shell matrix. More...

Detailed Description

Function Documentation

ConstrainMatrixDestroyOpPorQ()

MoFEMErrorCode ConstrainMatrixDestroyOpPorQ

(

Mat

Q

)

Destroy shell matrix Q.

  Mat Q; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
  CHKERR
MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ProjectionMatrixMultOpQ);
  CHKERR
MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpPorQ);

Definition at line 294 of file ConstrainMatrixCtx.cpp.

                                                   {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(Q, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  CHKERR ctx->destroyQorP();
  MoFEMFunctionReturn(0);
}

ConstrainMatrixDestroyOpQTKQ()

MoFEMErrorCode ConstrainMatrixDestroyOpQTKQ

(

Mat

QTKQ

)

Destroy shell matrix.

  Mat CTC_QTKQ; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
  CHKERR
MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpCTC_QTKQ);
  CHKERR
MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))mat_destroy_QTKQ);

Definition at line 302 of file ConstrainMatrixCtx.cpp.

                                                      {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(QTKQ, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  CHKERR ctx->destroyQTKQ();
  MoFEMFunctionReturn(0);
}

ConstrainMatrixMultOpCTC_QTKQ()

MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ	(	Mat	CTC_QTKQ,
		Vec	x,
		Vec	f
	)

Multiplication operator for RT = (CCT)^-TC.

  Mat CTC_QTKQ; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR
MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&CTC_QTKQ);
  CHKERR
MatShellSetOperation(CTC_QTKQ,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpCTC_QTKQ);
  CHKERR
MatShellSetOperation(CTC_QTKQ,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpQTKQ);

Definition at line 265 of file ConstrainMatrixCtx.cpp.

                                                                         {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(CTC_QTKQ, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  PetscLogEventBegin(ctx->MOFEM_EVENT_projCTC_QTKQ, 0, 0, 0, 0);
  Mat Q;
  int M, N, m, n;
  CHKERR MatGetSize(ctx->K, &M, &N);
  CHKERR MatGetLocalSize(ctx->K, &m, &n);
  CHKERR MatCreateShell(ctx->mField.get_comm(), m, n, M, N, ctx, &Q);
  CHKERR MatShellSetOperation(Q, MATOP_MULT,
                              (void (*)(void))ProjectionMatrixMultOpQ);
  CHKERR ctx->initializeQTKQ();
  CHKERR MatMult(Q, x, ctx->Qx);
  CHKERR MatMult(ctx->K, ctx->Qx, ctx->KQx);
  CHKERR MatMult(Q, ctx->KQx, f);
  CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
                         SCATTER_FORWARD);
  CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(ctx->CTC, ctx->X, ctx->CTCx);
  CHKERR VecScatterBegin(ctx->sCatter, ctx->CTCx, f, ADD_VALUES,
                         SCATTER_REVERSE);
  CHKERR VecScatterEnd(ctx->sCatter, ctx->CTCx, f, ADD_VALUES, SCATTER_REVERSE);
  CHKERR MatDestroy(&Q);
  PetscLogEventEnd(ctx->MOFEM_EVENT_projCTC_QTKQ, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}

ConstrainMatrixMultOpP()

MoFEMErrorCode ConstrainMatrixMultOpP	(	Mat	P,
		Vec	x,
		Vec	f
	)

Multiplication operator for P = CT(CCT)^-1C.

  Mat P; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&P);
  CHKERR
MatShellSetOperation(P,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpP);

Definition at line 201 of file ConstrainMatrixCtx.cpp.

                                                           {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(P, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  PetscLogEventBegin(ctx->MOFEM_EVENT_projP, 0, 0, 0, 0);
  CHKERR ctx->initializeQorP(x);
  CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
                         SCATTER_FORWARD);
  CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
  CHKERR KSPSolve(ctx->kSP, ctx->Cx, ctx->CCTm1_Cx);
  CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
  CHKERR VecZeroEntries(f);
  CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
                         SCATTER_REVERSE);
  CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
                       SCATTER_REVERSE);
  PetscLogEventEnd(ctx->MOFEM_EVENT_projP, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}

ConstrainMatrixMultOpR()

MoFEMErrorCode ConstrainMatrixMultOpR	(	Mat	R,
		Vec	x,
		Vec	f
	)

Multiplication operator for R = CT(CCT)^-1.

  Mat R; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&R);
  CHKERR
MatShellSetOperation(R,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpR);

Definition at line 225 of file ConstrainMatrixCtx.cpp.

                                                           {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(R, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  PetscLogEventBegin(ctx->MOFEM_EVENT_projR, 0, 0, 0, 0);
  if (ctx->initQorP)
    SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
            "you have to call first initQorP or use Q matrix");
  CHKERR KSPSolve(ctx->kSP, x, ctx->CCTm1_Cx);
  CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
  CHKERR VecZeroEntries(f);
  CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
                         SCATTER_REVERSE);
  CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
                       SCATTER_REVERSE);
  PetscLogEventEnd(ctx->MOFEM_EVENT_projR, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}

ConstrainMatrixMultOpRT()

MoFEMErrorCode ConstrainMatrixMultOpRT	(	Mat	RT,
		Vec	x,
		Vec	f
	)

Multiplication operator for RT = (CCT)^-TC.

  Mat RT; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&RT);
  CHKERR
MatShellSetOperation(RT,MATOP_MULT,(void(*)(void))ConstrainMatrixMultOpRT);

Definition at line 247 of file ConstrainMatrixCtx.cpp.

                                                             {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(RT, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  PetscLogEventBegin(ctx->MOFEM_EVENT_projRT, 0, 0, 0, 0);
  if (ctx->initQorP)
    SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
            "you have to call first initQorP or use Q matrix");
  CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
                         SCATTER_FORWARD);
  CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
  CHKERR KSPSolve(ctx->kSP, ctx->Cx, f);
  PetscLogEventEnd(ctx->MOFEM_EVENT_projRT, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}

ProjectionMatrixMultOpQ()

MoFEMErrorCode ProjectionMatrixMultOpQ	(	Mat	Q,
		Vec	x,
		Vec	f
	)

Multiplication operator for Q = I-CTC(CCT)^-1C.

  Mat Q; //for problem
  ConstrainMatrixCtx
projection_matrix_ctx(m_field,problem_name,contrains_problem_name);
  CHKERR MatCreateShell(PETSC_COMM_WORLD,m,m,M,M,&projection_matrix_ctx,&Q);
  CHKERR
MatShellSetOperation(Q,MATOP_MULT,(void(*)(void))ProjectionMatrixMultOpQ);
  CHKERR
MatShellSetOperation(Q,MATOP_DESTROY,(void(*)(void))ConstrainMatrixDestroyOpPorQ);

Definition at line 167 of file ConstrainMatrixCtx.cpp.

                                                            {
  MoFEMFunctionBegin;
  void *void_ctx;
  CHKERR MatShellGetContext(Q, &void_ctx);
  ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
  PetscLogEventBegin(ctx->MOFEM_EVENT_projQ, 0, 0, 0, 0);
  CHKERR ctx->initializeQorP(x);
  CHKERR VecCopy(x, f);
  CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
                         SCATTER_FORWARD);
  CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
  if (debug) {
    // CHKERR VecView(ctx->X,PETSC_VIEWER_STDOUT_WORLD);
    CHKERR VecScatterBegin(ctx->sCatter, ctx->X, f, INSERT_VALUES,
                           SCATTER_REVERSE);
    CHKERR VecScatterEnd(ctx->sCatter, ctx->X, f, INSERT_VALUES,
                         SCATTER_REVERSE);
    PetscBool flg;
    CHKERR VecEqual(x, f, &flg);
    if (flg == PETSC_FALSE)
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "scatter is not working");
  }
  CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
  CHKERR KSPSolve(ctx->kSP, ctx->Cx, ctx->CCTm1_Cx);
  CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
  CHKERR VecScale(ctx->CT_CCTm1_Cx, -1);
  CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, ADD_VALUES,
                         SCATTER_REVERSE);
  CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, ADD_VALUES,
                       SCATTER_REVERSE);
  PetscLogEventEnd(ctx->MOFEM_EVENT_projQ, 0, 0, 0, 0);
  MoFEMFunctionReturn(0);
}