C++ NumericVector::add_vector_blocked方法代码示例

//.........这里部分代码省略.........
		  else if (2 == ivar + jvar ) kvar = dim+2; // xz
		  else if (3 == ivar + jvar ) kvar = dim+1; // yz
		  Res[ivar] += - SolVAR[jvar]*GradSolVAR[ivar][jvar] 
			       + IRe * ( NablaSolVAR[ivar][jvar] + NablaSolVAR[jvar][kvar] );
		}
	      }
	    
	      adept::adouble div_vel=0.;
	      for(int ivar=0; ivar<dim; ivar++) {
		div_vel +=GradSolVAR[ivar][ivar];
	      }
	    	    
	      //BEGIN redidual momentum block  
	      for (unsigned i=0; i<nve2; i++){
		for(unsigned ivar=0; ivar<dim; ivar++) {
				
		  adept::adouble Advection = 0.;
		  adept::adouble Laplacian = 0.;
		  adept::adouble phiSupg=0.;
		  for(unsigned jvar=0; jvar<dim; jvar++) {
		    unsigned kvar;
		    if(ivar == jvar) kvar = jvar;
		    else if (1 == ivar + jvar ) kvar = dim;   // xy
		    else if (2 == ivar + jvar ) kvar = dim+2; // xz
		    else if (3 == ivar + jvar ) kvar = dim+1; // yz
		
		    Advection += SolVAR[jvar]*GradSolVAR[ivar][jvar]*phi[i];
		    Laplacian += IRe*gradphi[i*dim+jvar]*(GradSolVAR[ivar][jvar]+GradSolVAR[jvar][ivar]);
		    phiSupg   += ( SolVAR[jvar]*gradphi[i*dim+jvar] ) * tauSupg;  
		    
		    aRhs[indexVAR[ivar]][i] += 	 Res[ivar] * (-IRe * nablaphi[i*nabla_dim+jvar])* tauSupg * Weight; 
		    aRhs[indexVAR[jvar]][i] += 	 Res[ivar] * (-IRe * nablaphi[i*nabla_dim+kvar])* tauSupg * Weight; 		  
		  }
		  aRhs[indexVAR[ivar]][i]+= ( - Advection - Laplacian 
					      + ( SolVAR[dim] - deltaSupg * div_vel) * gradphi[i*dim+ivar]
					      + Res[ivar] * phiSupg ) * Weight;      
		}
	      } 
	      //END redidual momentum block     
	    
	      //BEGIN continuity block 
	      for (unsigned i=0; i<nve1; i++) {
		adept::adouble MinusGradPhi1DotRes = 0.;
		for(int ivar=0;ivar<dim;ivar++){
		  MinusGradPhi1DotRes += -gradphi1[i*dim+ivar] * Res[ivar] * tauSupg; 
		}
		aRhs[indexVAR[dim]][i] += ( - (-div_vel) * phi1[i] + MinusGradPhi1DotRes ) * Weight;
	      }
	      //END continuity block ===========================
	    }	   
	    //END FLUID ASSEMBLY ============
	  }
	}
      }
	
      //BEGIN local to global assembly 	
      //copy adouble aRhs into double Rhs
      for (unsigned i=0;i<dim;i++) {
	for(int j=0; j<nve2; j++) {
	  Rhs[indexVAR[i]][j] = aRhs[indexVAR[i]][j].value();
	}
      }
      for (unsigned j=0;j<nve1;j++) {
	Rhs[indexVAR[dim]][j] = aRhs[indexVAR[dim]][j].value();
      }	
      for(int i=0; i<dim+1; i++) {
	myRES->add_vector_blocked(Rhs[indexVAR[i]],dofsVAR[i]);
      }     
                  
      //Store equations
      for(int i=0; i<dim; i++) {  
	adeptStack.dependent(&aRhs[indexVAR[i]][0], nve2);
	adeptStack.independent(&Soli[indexVAR[i]][0], nve2); 
      }
      adeptStack.dependent(&aRhs[indexVAR[dim]][0], nve1);
      adeptStack.independent(&Soli[indexVAR[dim]][0], nve1);   
      adeptStack.jacobian(&Jac[0]);	
      unsigned nveAll=(dim*nve2+nve1);
      for (int inode=0;inode<nveAll;inode++){
	for (int jnode=0;jnode<nveAll;jnode++){
	   KKloc[inode*nveAll+jnode]=-Jac[jnode*nveAll+inode];
	}
      }
      myKK->add_matrix_blocked(KKloc,dofsAll,dofsAll);
      adeptStack.clear_independents();
      adeptStack.clear_dependents();
       
      //END local to global assembly
   
    } //end list of elements loop
    
    myKK->close();
    myRES->close();
           
    // *************************************
    end_time=clock();
    AssemblyTime+=(end_time-start_time);
    // ***************** END ASSEMBLY RESIDUAL + MATRIX *******************

  }  

//.........这里部分代码省略.........
      solu[i]          = (*sol->_Sol[soluIndex])(solDof);      // global extraction and local storage for the solution
      solv[i]          = (*sol->_Sol[solvIndex])(solDof);      // global extraction and local storage for the solution
      sysDof[i]         = pdeSys->GetSystemDof(soluIndex, soluPdeIndex, i, iel);    // global to global mapping between solution node and pdeSys dof
      sysDof[nDofs + i] = pdeSys->GetSystemDof(solvIndex, solvPdeIndex, i, iel);    // global to global mapping between solution node and pdeSys dof
    }

    // local storage of coordinates
    for (unsigned i = 0; i < nDofs2; i++) {
      unsigned xDof  = msh->GetSolutionDof(i, iel, xType); // global to global mapping between coordinates node and coordinate dof

      for (unsigned jdim = 0; jdim < dim; jdim++) {
        x[jdim][i] = (*msh->_topology->_Sol[jdim])(xDof);  // global extraction and local storage for the element coordinates
      }
    }

    // start a new recording of all the operations involving adept::adouble variables
    s.new_recording();

    // *** Gauss point loop ***
    for (unsigned ig = 0; ig < msh->_finiteElement[ielGeom][soluType]->GetGaussPointNumber(); ig++) {
      // *** get gauss point weight, test function and test function partial derivatives ***
      msh->_finiteElement[ielGeom][soluType]->Jacobian(x, ig, weight, phi, phi_x, phi_xx);

      // evaluate the solution, the solution derivatives and the coordinates in the gauss point
      adept::adouble soluGauss = 0;
      vector < adept::adouble > soluGauss_x(dim, 0.);

      adept::adouble solvGauss = 0;
      vector < adept::adouble > solvGauss_x(dim, 0.);

      vector < double > xGauss(dim, 0.);

      for (unsigned i = 0; i < nDofs; i++) {
        soluGauss += phi[i] * solu[i];
        solvGauss += phi[i] * solv[i];

        for (unsigned jdim = 0; jdim < dim; jdim++) {
          soluGauss_x[jdim] += phi_x[i * dim + jdim] * solu[i];
          solvGauss_x[jdim] += phi_x[i * dim + jdim] * solv[i];
          xGauss[jdim] += x[jdim][i] * phi[i];
        }
      }

      // *** phi_i loop ***
      for (unsigned i = 0; i < nDofs; i++) {

        adept::adouble Laplace_u = 0.;
        adept::adouble Laplace_v = 0.;

        for (unsigned jdim = 0; jdim < dim; jdim++) {
          Laplace_u   +=  - phi_x[i * dim + jdim] * soluGauss_x[jdim];
          Laplace_v   +=  - phi_x[i * dim + jdim] * solvGauss_x[jdim];
        }

        double exactSolValue = GetExactSolutionValue(xGauss);

        double pi = acos(-1.);
        aResv[i] += (solvGauss * phi[i] -  Laplace_u) * weight;
        aResu[i] += (4.*pi * pi * pi * pi * exactSolValue * phi[i] -  Laplace_v) * weight;
      } // end phi_i loop
    } // end gauss point loop


    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store
    Res.resize(2 * nDofs);

    for (int i = 0; i < nDofs; i++) {
      Res[i]         = -aResu[i].value();
      Res[nDofs + i] = -aResv[i].value();
    }

    RES->add_vector_blocked(Res, sysDof);

    Jac.resize( 4 * nDofs * nDofs );

    // define the dependent variables
    s.dependent(&aResu[0], nDofs);
    s.dependent(&aResv[0], nDofs);

    // define the independent variables
    s.independent(&solu[0], nDofs);
    s.independent(&solv[0], nDofs);

    // get the jacobian matrix (ordered by column)
    s.jacobian(&Jac[0], true);

    KK->add_matrix_blocked(Jac, sysDof, sysDof);

    s.clear_independents();
    s.clear_dependents();

  } //end element loop for each process

  RES->close();
  KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
      msh->_finiteElement[ielGeom][solUType]->Jacobian (crdX, ig, weight, phi, phi_x, phi_xx);
      msh->_finiteElement[ielGeom][solVType]->Jacobian (crdX, ig, weight, phiV, phiV_x, phiV_xx);

      adept::adouble solUig = 0; // solution U in the gauss point
      vector < adept::adouble > gradSolUig (dim, 0.); // gradient of solution U in the gauss point

      for (unsigned i = 0; i < nDofsU; i++) {
        solUig += phi[i] * solU[i];

        for (unsigned j = 0; j < dim; j++) {
          gradSolUig[j] += phi_x[i * dim + j] * solU[i];
        }
      }

      adept::adouble solVig = 0; // solution V in the gauss point
      vector < adept::adouble > gradSolVig (dim, 0.); // gradient of solution U in the gauss point

      for (unsigned i = 0; i < nDofsV; i++) {
        solVig += phiV[i] * solV[i];

        for (unsigned j = 0; j < dim; j++) {
          gradSolVig[j] += phiV_x[i * dim + j] * solV[i];
        }
      }

      double nu = 1.;

      // *** phiU_i loop ***
      for (unsigned i = 0; i < nDofsU; i++) {
        adept::adouble LaplaceU = 0.;

        for (unsigned j = 0; j < dim; j++) {
          LaplaceU -=  nu * phi_x[i * dim + j] * gradSolUig[j];
        }

        aResU[i] += (phi[i] - LaplaceU) * weight;

      } // end phiU_i loop

      for (unsigned i = 0; i < nDofsV; i++) {
        adept::adouble LaplaceV = 0.;

        for (unsigned j = 0; j < dim; j++) {
          LaplaceV -=  nu * phiV_x[i * dim + j] * gradSolVig[j];
        }

        aResV[i] += (phiV[i] * solUig - LaplaceV) * weight;

      } // end phiV_i loop

    } // end gauss point loop




    // } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector



    //copy the value of the adept::adoube aRes in double Res and store them in RES
    Res.resize (nDofsU + nDofsV);   //resize

    for (int i = 0; i < nDofsU; i++) {
      Res[i] = -aResU[i].value();
    }

    for (int i = 0; i < nDofsV; i++) {
      Res[i + nDofsU] = -aResV[i].value();
    }

    RES->add_vector_blocked (Res, sysDof);

    //Extarct and store the Jacobian

    Jac.resize ( (nDofsU + nDofsV) * (nDofsU + nDofsV));
    // define the dependent variables
    s.dependent (&aResU[0], nDofsU);
    s.dependent (&aResV[0], nDofsV);

    // define the independent variables
    s.independent (&solU[0], nDofsU);
    s.independent (&solV[0], nDofsV);
    // get the and store jacobian matrix (row-major)
    s.jacobian (&Jac[0] , true);
    KK->add_matrix_blocked (Jac, sysDof, sysDof);

    s.clear_independents();
    s.clear_dependents();

  } //end element loop for each process

  RES->close();

  KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
          for (unsigned  k = 0; k < dim; k++) {
            solV_gss[k] += phiV[i] * solV[k][i];
          }

          for (unsigned j = 0; j < dim; j++) {
            for (unsigned  k = 0; k < dim; k++) {
              gradSolV_gss[k][j] += phiV_x[i * dim + j] * solV[k][i];
            }
          }
        }

        adept::adouble solP_gss = 0;

        for (unsigned i = 0; i < nDofsP; i++) {
          solP_gss += phiP[i] * solP[i];
        }

        double nu = 1.;

        // *** phiV_i loop ***
        for (unsigned i = 0; i < nDofsV; i++) {
          vector < adept::adouble > NSV(dim, 0.);

          for (unsigned j = 0; j < dim; j++) {
            for (unsigned  k = 0; k < dim; k++) {
              NSV[k]   +=  nu * phiV_x[i * dim + j] * (gradSolV_gss[k][j] + gradSolV_gss[j][k]);
              NSV[k]   +=  phiV[i] * (solV_gss[j] * gradSolV_gss[k][j]);
            }
          }

          for (unsigned  k = 0; k < dim; k++) {
            NSV[k] += -solP_gss * phiV_x[i * dim + k];
          }

          for (unsigned  k = 0; k < dim; k++) {
            aResV[k][i] += - NSV[k] * weight;
          }
        } // end phiV_i loop

        // *** phiP_i loop ***
        for (unsigned i = 0; i < nDofsP; i++) {
          for (int k = 0; k < dim; k++) {
            aResP[i] += - (gradSolV_gss[k][k]) * phiP[i]  * weight;
          }
        } // end phiP_i loop

      } // end gauss point loop
    } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store them in RES
    Res.resize(nDofsVP);    //resize

    for (int i = 0; i < nDofsV; i++) {
      for (unsigned  k = 0; k < dim; k++) {
        Res[ i +  k * nDofsV ] = -aResV[k][i].value();
      }
    }

    for (int i = 0; i < nDofsP; i++) {
      Res[ i + dim * nDofsV ] = -aResP[i].value();
    }

    RES->add_vector_blocked(Res, KKDof);

    //Extarct and store the Jacobian
    if (assembleMatrix) {
      Jac.resize(nDofsVP * nDofsVP);
      // define the dependent variables

      for (unsigned  k = 0; k < dim; k++) {
        s.dependent(&aResV[k][0], nDofsV);
      }

      s.dependent(&aResP[0], nDofsP);

      // define the independent variables
      for (unsigned  k = 0; k < dim; k++) {
        s.independent(&solV[k][0], nDofsV);
      }

      s.independent(&solP[0], nDofsP);

      // get the and store jacobian matrix (row-major)
      s.jacobian(&Jac[0] , true);
      KK->add_matrix_blocked(Jac, KKDof, KKDof);

      s.clear_independents();
      s.clear_dependents();
    }
  } //end element loop for each process

  RES->close();

  if (assembleMatrix) KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........

        }
      }

      //END phiT_i loop

      //BEGIN phiV_i loop: Momentum balance
      for(unsigned i = 0; i < nDofsV; i++) {

        for(unsigned k = 0; k < dim; k++) {
          unsigned irow = nDofsT + k * nDofsV + i;

          for(unsigned l = 0; l < dim; l++) {
            Res[irow] +=  -sqrt(Pr / Ra) * phiV_x[i * dim + l] * (gradSolV_gss[k][l] + gradSolV_gss[l][k]) * weight;
            Res[irow] +=  -phiV[i] * solV_gss[l] * gradSolV_gss[k][l] * weight;
          }

          Res[irow] += solP_gss * phiV_x[i * dim + k] * weight;

          if(k == 1) {
            Res[irow] += beta * solT_gss * phiV[i] * weight;
          }

          if(assembleMatrix) {
            unsigned irowMat = nDofsTVP * irow;

            for(unsigned l = 0; l < dim; l++) {
              for(unsigned j = 0; j < nDofsV; j++) {
                unsigned jcol1 = (nDofsT + k * nDofsV + j);
                unsigned jcol2 = (nDofsT + l * nDofsV + j);
                Jac[ irowMat + jcol1] += sqrt(Pr / Ra) * phiV_x[i * dim + l] * phiV_x[j * dim + l] * weight;
                Jac[ irowMat + jcol2] += sqrt(Pr / Ra) * phiV_x[i * dim + l] * phiV_x[j * dim + k] * weight;
                Jac[ irowMat + jcol1] += phiV[i] * solV_gss[l] * phiV_x[j * dim + l] * weight;
                Jac[ irowMat + jcol2] += phiV[i] * phiV[j] * gradSolV_gss[k][l] * weight;
              }
            }

            for(unsigned j = 0; j < nDofsP; j++) {
              unsigned jcol = (nDofsT + dim * nDofsV) + j;
              Jac[ irowMat + jcol] += - phiV_x[i * dim + k] * phiP[j] * weight;
            }

            if(k == 1) {
              for(unsigned j = 0; j < nDofsT; j++) {
                Jac[ irowMat + j ] += -beta * phiV[i] * phiT[j] * weight;
              }
            }
          }

        }
      }

      //END phiV_i loop

      //BEGIN phiP_i loop: mass balance
      for(unsigned i = 0; i < nDofsP; i++) {
        unsigned irow = nDofsT + dim * nDofsV + i;

        for(int k = 0; k < dim; k++) {
          Res[irow] += (gradSolV_gss[k][k]) * phiP[i]  * weight;

          if(assembleMatrix) {
            unsigned irowMat = nDofsTVP * irow;

            for(unsigned j = 0; j < nDofsV; j++) {
              unsigned jcol = (nDofsT + k * nDofsV + j);
              Jac[ irowMat + jcol ] += - phiP[i] * phiV_x[j * dim + k] * weight;
            }
          }

        }
      }

      //END phiP_i loop

    }

    //END Gauss point loop

    //BEGIN local to global Matrix/Vector assembly
    RES->add_vector_blocked(Res, sysDof);

    if(assembleMatrix) {
      KK->add_matrix_blocked(Jac, sysDof, sysDof);
    }

    //END local to global Matrix/Vector assembly

  }

  //END element loop

  RES->close();

  if(assembleMatrix) {
    KK->close();
  }

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
        for (unsigned j = 0; j < nDofv; j++) {
          double sign = ( j == 0) ? 1. : -1;
          aResHT[k][0] += sign * solv[k][j] * solHT[k][0] / dx; 
	}
	if(k<NLayers-1){
	  aResHT[k][0] += w[k+1] * 0.5 * (solHT[k][0]/solh[k][0] + solHT[k+1][0]/solh[k+1][0]);
	  //aResHT[k][0] += w[k+1] * 0.5 * (solT[k][0] + solT[k+1][0]);
	}
	if( k > 0){
	  aResHT[k][0] -= w[k] * 0.5 * (solHT[k-1][0]/solh[k-1][0] + solHT[k][0]/solh[k][0] );
	  //aResHT[k][0] -= w[k] * 0.5 * (solT[k-1][0] + solT[k][0]);
	}
      }
    } 

    vector< double > Res(NLayers * nDofs); // local redidual vector

    unsigned counter = 0;
    for(unsigned k = 0; k < NLayers; k++) {
      for(int i = 0; i < nDofh; i++) {
        Res[counter] =  aResh[k][i].value();
        counter++;
      }
      for(int i = 0; i < nDofv; i++) {
        Res[counter] =  aResv[k][i].value();
        counter++;
      }
      for(int i = 0; i < nDofHT; i++) {
        Res[counter] =  aResHT[k][i].value();
        counter++;
      }
    }

    RES->add_vector_blocked(Res, l2GMap);


    for(unsigned k = 0; k < NLayers; k++) {
      // define the dependent variables
      s.dependent(&aResh[k][0], nDofh);
      s.dependent(&aResv[k][0], nDofv);
      s.dependent(&aResHT[k][0], nDofHT);

      // define the independent variables
      s.independent(&solh[k][0], nDofh);
      s.independent(&solv[k][0], nDofv);
      s.independent(&solHT[k][0], nDofHT);
    }

    // get the jacobian matrix (ordered by row major )
    vector < double > Jac(NLayers * nDofs * NLayers * nDofs);
    s.jacobian(&Jac[0], true);

    //store K in the global matrix KK
    KK->add_matrix_blocked(Jac, l2GMap, l2GMap);

    s.clear_independents();
    s.clear_dependents();

  }

  RES->close();
  KK->close();

//   PetscViewer    viewer;
//   PetscViewerDrawOpen(PETSC_COMM_WORLD,NULL,NULL,0,0,900,900,&viewer);
//   PetscObjectSetName((PetscObject)viewer,"FSI matrix");

//.........这里部分代码省略.........
            const unsigned felt = mymsh->GetElementFaceType(iel, jface);

            for (unsigned i = 0; i < nve; i++) {
              unsigned ilocal = mymsh->GetLocalFaceVertexIndex(iel, jface, i);
              unsigned inode_coord_metis = mymsh->GetSolutionDof(ilocal, iel, 2);

              for (unsigned ivar = 0; ivar < dim; ivar++) {
                coordinates[ivar][i] = (*mymsh->_topology->_Sol[ivar])(inode_coord_metis);
              }
            }

            if (felt != 6) {
              for (unsigned igs = 0; igs < mymsh->_finiteElement[felt][order_ind]->GetGaussPointNumber(); igs++) {
                mymsh->_finiteElement[felt][order_ind]->JacobianSur(coordinates, igs, weight, phi, gradphi, normal);

                xyzt.assign(4, 0.);

                for (unsigned i = 0; i < nve; i++) {
                  for (unsigned ivar = 0; ivar < dim; ivar++) {
                    xyzt[ivar] += coordinates[ivar][i] * phi[i];
                  }
                }

                // *** phi_i loop ***
                for (unsigned i = 0; i < nve; i++) {
                  double surfterm_g = (*bdcfunc)(&xyzt[0]);
                  double bdintegral = phi[i] * surfterm_g * weight;
                  unsigned int ilocalnode = mymsh->GetLocalFaceVertexIndex(iel, jface, i);
                  F[ilocalnode] += bdintegral;
                }
              }
            }
            else { // 1D : the side elems are points and does not still exist the point elem
              // in 1D it is only one point
              xyzt[0] = coordinates[0][0];
              xyzt[1] = 0.;
              xyzt[2] = 0.;
              xyzt[3] = 0.;

              double bdintegral = (*bdcfunc)(&xyzt[0]);
              unsigned int ilocalnode = mymsh->GetLocalFaceVertexIndex(iel, jface, 0);
              F[ilocalnode] += bdintegral;
            }
          }
        }
      }
    }
    else {
      double tau = 0.;
      std::vector< double > xx(dim, 0.);
      vector < double > normal(dim, 0);

      // loop on faces
      for (unsigned jface = 0; jface < mymsh->GetElementFaceNumber(iel); jface++) {
        // look for boundary faces
        if (myel->GetFaceElementIndex(iel, jface) < 0) {
          unsigned int faceIndex =  myel->GetBoundaryIndex(iel, jface);

          if (!ml_sol->GetBdcFunction()(xx, "Sol", tau, faceIndex, 0.) && tau != 0.) {
            unsigned nve = mymsh->GetElementFaceDofNumber(iel, jface, order_ind);
            const unsigned felt = mymsh->GetElementFaceType(iel, jface);

            for (unsigned i = 0; i < nve; i++) {
              unsigned int ilocal = mymsh->GetLocalFaceVertexIndex(iel, jface, i);
              unsigned inode = mymsh->GetSolutionDof(ilocal, iel, 2);

              for (unsigned idim = 0; idim < dim; idim++) {
                coordinates[idim][i] = (*mymsh->_topology->_Sol[idim])(inode);
              }
            }

            for (unsigned igs = 0; igs < mymsh->_finiteElement[felt][order_ind]->GetGaussPointNumber(); igs++) {
              mymsh->_finiteElement[felt][order_ind]->JacobianSur(coordinates, igs, weight, phi, gradphi, normal);

              // *** phi_i loop ***
              for (unsigned i = 0; i < nve; i++) {
                double value = phi[i] * tau * weight;
                unsigned int ilocal = mymsh->GetLocalFaceVertexIndex(iel, jface, i);
                F[ilocal]   += value;
              }
            }
          }
        }
      }
    }

    //--------------------------------------------------------------------------------------------------------
    //Sum the local matrices/vectors into the global Matrix/vector

    myRES->add_vector_blocked(F, KK_dof);

    myKK->add_matrix_blocked(B, KK_dof, KK_dof);
  } //end list of elements loop for each subdomain

  myRES->close();
  myKK->close();

  // ***************** END ASSEMBLY *******************

}

//.........这里部分代码省略.........
    }

    Res.resize(nDofs);    //resize
    std::fill(Res.begin(), Res.end(), 0);    //set Res to zero

    if (assembleMatrix) {
      J.resize(nDofs * nDofs);    //resize
      std::fill(J.begin(), J.end(), 0);    //set K to zero
    }

    // local storage of global mapping and solution
    for (unsigned i = 0; i < nDofs; i++) {
      unsigned iNode = el->GetMeshDof(kel, i, soluType);    // local to global solution node
      unsigned solDof = msh->GetMetisDof(iNode, soluType);    // global to global mapping between solution node and solution dof
      solu[i] = (*sol->_Sol[soluIndex])(solDof);      // global extraction and local storage for the solution
      KKDof[i] = pdeSys->GetKKDof(soluIndex, soluPdeIndex, iNode);    // global to global mapping between solution node and pdeSys dof
    }

    // local storage of coordinates
    for (unsigned i = 0; i < nDofs2; i++) {
      unsigned iNode = el->GetMeshDof(kel, i, xType);    // local to global coordinates node
      unsigned xDof  = msh->GetMetisDof(iNode, xType);    // global to global mapping between coordinates node and coordinate dof

      for (unsigned jdim = 0; jdim < dim; jdim++) {
        x[jdim][i] = (*msh->_coordinate->_Sol[jdim])(xDof);      // global extraction and local storage for the element coordinates
      }
    }

    if (level == levelMax || !el->GetRefinedElementIndex(kel)) {      // do not care about this if now (it is used for the AMR)
      // *** Gauss point loop ***
      for (unsigned ig = 0; ig < msh->_finiteElement[kelGeom][soluType]->GetGaussPointNumber(); ig++) {
        // *** get gauss point weight, test function and test function partial derivatives ***
        msh->_finiteElement[kelGeom][soluType]->Jacobian(x, ig, weight, phi, phi_x, phi_xx);

        // evaluate the solution, the solution derivatives and the coordinates in the gauss point
        double soluGauss = 0;
        vector < double > soluGauss_x(dim, 0.);
        vector < double > xGauss(dim, 0.);

        for (unsigned i = 0; i < nDofs; i++) {
          soluGauss += phi[i] * solu[i];

          for (unsigned jdim = 0; jdim < dim; jdim++) {
            soluGauss_x[jdim] += phi_x[i * dim + jdim] * solu[i];
            xGauss[jdim] += x[jdim][i] * phi[i];
          }
        }

        // *** phi_i loop ***
        for (unsigned i = 0; i < nDofs; i++) {

          double nonLinearTerm = 0.;
          double mLaplace = 0.;

          for (unsigned jdim = 0; jdim < dim; jdim++) {
            mLaplace   +=  phi_x[i * dim + jdim] * soluGauss_x[jdim];
            nonLinearTerm += soluGauss * soluGauss_x[jdim] * phi[i];
          }

          double exactSolValue = GetExactSolutionValue(xGauss);
          vector < double > exactSolGrad(dim);
          GetExactSolutionGradient(xGauss , exactSolGrad);
          double exactSolLaplace = GetExactSolutionLaplace(xGauss);


          double f = (- exactSolLaplace + exactSolValue * (exactSolGrad[0] + exactSolGrad[1])) * phi[i] ;
          Res[i] += (f - (mLaplace + nonLinearTerm)) * weight;

          if (assembleMatrix) {
            // *** phi_j loop ***
            for (unsigned j = 0; j < nDofs; j++) {
              mLaplace = 0.;
              nonLinearTerm = 0.;

              for (unsigned kdim = 0; kdim < dim; kdim++) {
                mLaplace += (phi_x[i * dim + kdim] * phi_x[j * dim + kdim]);
                nonLinearTerm +=   phi[i] * (phi[j] * soluGauss_x[kdim] + soluGauss * phi_x[j * dim + kdim]);
              }

              J[i * nDofs + j] += (mLaplace + nonLinearTerm) * weight;
            } // end phi_j loop
          } // endif assemble_matrix
        } // end phi_i loop
      } // end gauss point loop
    } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    RES->add_vector_blocked(Res, KKDof);

    if (assembleMatrix) KK->add_matrix_blocked(J, KKDof, KKDof);
  } //end element loop for each process

  RES->close();

  if (assembleMatrix) KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........

            h[0][0] = h[0][1] = h[1][0] = h[1][1] = 0.;

            for (int k = 0; k < 3; k++) {
                for (int u = 0; u < dim; u++) {
                    for (int v = 0; v < dim; v++) {
                        h[u][v] += solRGauss_xx[k][u][v] * N[k];

                    }
                }
            }

            //adept::adouble K = cos(sol_x[0])/(a+cos(sol_x[0]));//(h[0][0]*h[1][1]-h[0][1]*h[1][0])/detg;

            adept::adouble K = (h[0][0] * h[1][1] - h[0][1] * h[1][0]) / detg;

            adept::adouble H_exact = 0.5 * (1. + cos(sol_x[0]) / (a + cos(sol_x[0])));

            // *** phi_i loop ***
            for (unsigned i = 0; i < nDofs; i++) {

                for (int k = 0; k < 3; k++) {
                    for (int u = 0; u < dim; u++) {
                        adept::adouble AgIgradRgradPhi = 0;

                        for (int v = 0; v < dim; v++) {
                            AgIgradRgradPhi += A * gI[u][v].value() * solRGauss_x[k][v];
                        }

                        aResR[k][i] += AgIgradRgradPhi * phi_x[i * dim + u] * weight;
                    }

                    aResR[k][i] += 2.* A * solHGauss.value() * N[k] * phi[i] * weight;
                }


                for (int u = 0; u < dim; u++) {
                    adept::adouble AgIgradHgradPhi = 0;

                    for (int v = 0; v < dim; v++) {
                        AgIgradHgradPhi += A * gI[u][v].value() * solHGauss_x[v];
                    }

                    aResH[i] -= AgIgradHgradPhi * phi_x[i * dim + u] * weight;
                }

                aResH[i] += A * (-0 * (solRGauss[0] - solRGaussOld[0]) * N[0].value()
                                 - 0 * (solRGauss[1] - solRGaussOld[1]) * N[1].value()
                                 - 0 * (solRGauss[2] - solRGaussOld[2]) * N[2].value()
                                 + 2. * solHGauss * (solHGauss * solHGauss  - K.value())) * phi[i] * weight;

            } // end phi_i loop
        } // end gauss point loop

        //--------------------------------------------------------------------------------------------------------
        // Add the local Matrix/Vector into the global Matrix/Vector

        //copy the value of the adept::adoube aRes in double Res and store
        for (int i = 0; i < nDofs; i++) {
            for (int k = 0; k < 3; k++) {
                Res[ k * nDofs + i] = -aResR[k][i].value();
            }

            Res[ 3 * nDofs + i] = -aResH[i].value();
        }

        RES->add_vector_blocked(Res, sysDof);


        Jac.resize((4 * nDofs) * (4 * nDofs));

        // define the dependent variables
        for (int k = 0; k < 3; k++) {
            s.dependent(&aResR[k][0], nDofs);
        }

        s.dependent(&aResH[0], nDofs);

        // define the independent variables
        for (int k = 0; k < 3; k++) {
            s.independent(&solR[k][0], nDofs);
        }

        s.independent(&solH[0], nDofs);
        // get the jacobian matrix (ordered by row)
        s.jacobian(&Jac[0], true);

        KK->add_matrix_blocked(Jac, sysDof, sysDof);

        s.clear_independents();
        s.clear_dependents();

    } //end element loop for each process

    RES->close();

    KK->close();

    // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........

    for (int i = 0; i < dim; i++) {
      x[i].resize(nDofx);
    }

    aRes.resize(nDofu);    //resize
    std::fill(aRes.begin(), aRes.end(), 0);    //set aRes to zero

    // local storage of global mapping and solution
    for (unsigned i = 0; i < nDofu; i++) {
      unsigned solDof = msh->GetSolutionDof(i, iel, soluType);    // global to global mapping between solution node and solution dof
      solu[i] = (*sol->_Sol[soluIndex])(solDof);      // global extraction and local storage for the solution
      l2GMap[i] = pdeSys->GetSystemDof(soluIndex, soluPdeIndex, i, iel);    // global to global mapping between solution node and pdeSys dof
    }

    // local storage of coordinates
    for (unsigned i = 0; i < nDofx; i++) {
      unsigned xDof  = msh->GetSolutionDof(i, iel, xType);    // global to global mapping between coordinates node and coordinate dof

      for (unsigned jdim = 0; jdim < dim; jdim++) {
        x[jdim][i] = (*msh->_topology->_Sol[jdim])(xDof);      // global extraction and local storage for the element coordinates
      }
    }


    // start a new recording of all the operations involving adept::adouble variables
    s.new_recording();

    // *** Gauss point loop ***
    for (unsigned ig = 0; ig < msh->_finiteElement[ielGeom][soluType]->GetGaussPointNumber(); ig++) {
      // *** get gauss point weight, test function and test function partial derivatives ***
      msh->_finiteElement[ielGeom][soluType]->Jacobian(x, ig, weight, phi, phi_x, phi_xx);

      // evaluate the solution, the solution derivatives and the coordinates in the gauss point
      adept::adouble solu_gss = 0;
      vector < adept::adouble > gradSolu_gss(dim, 0.);
      vector < double > x_gss(dim, 0.);

      for (unsigned i = 0; i < nDofu; i++) {
        solu_gss += phi[i] * solu[i];

        for (unsigned jdim = 0; jdim < dim; jdim++) {
          gradSolu_gss[jdim] += phi_x[i * dim + jdim] * solu[i];
          x_gss[jdim] += x[jdim][i] * phi[i];
        }
      }

      // *** phi_i loop ***
      for (unsigned i = 0; i < nDofu; i++) {

        adept::adouble laplace = 0.;

        for (unsigned jdim = 0; jdim < dim; jdim++) {
          laplace   +=  phi_x[i * dim + jdim] * gradSolu_gss[jdim];
        }

        double srcTerm = - GetExactSolutionLaplace(x_gss);
        aRes[i] += (srcTerm * phi[i] - laplace) * weight;

      } // end phi_i loop
    } // end gauss point loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store
    Res.resize(nDofu);    //resize

    for (int i = 0; i < nDofu; i++) {
      Res[i] = - aRes[i].value();
    }

    RES->add_vector_blocked(Res, l2GMap);



    // define the dependent variables
    s.dependent(&aRes[0], nDofu);

    // define the independent variables
    s.independent(&solu[0], nDofu);

    // get the jacobian matrix (ordered by row major )
    Jac.resize(nDofu * nDofu);    //resize
    s.jacobian(&Jac[0], true);

    //store K in the global matrix KK
    KK->add_matrix_blocked(Jac, l2GMap, l2GMap);

    s.clear_independents();
    s.clear_dependents();

  } //end element loop for each process

  RES->close();

  KK->close();

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........


      // *** phiV_i loop ***
      for(unsigned i = 0; i < nDofsV; i++) {
        vector < adept::adouble > NSV(dim, 0.);
        vector < double > NSVold(dim, 0.);

        for(unsigned j = 0; j < dim; j++) {
          for(unsigned  k = 0; k < dim; k++) {
            NSV[k]   +=  sqrt(Pr / Ra) * phiV_x[i * dim + j] * (gradSolV_gss[k][j] + gradSolV_gss[j][k]);
            NSV[k]   +=  phiV[i] * (solV_gss[j] * gradSolV_gss[k][j]);

            NSVold[k]   +=  sqrt(Pr / Ra) * phiV_x[i * dim + j] * (gradSolVold_gss[k][j] + gradSolVold_gss[j][k]);
            NSVold[k]   +=  phiV[i] * (solVold_gss[j] * gradSolVold_gss[k][j]);

          }
        }

        for(unsigned  k = 0; k < dim; k++) {
          NSV[k] += -solP_gss * phiV_x[i * dim + k];
          NSVold[k] += -solPold_gss * phiV_x[i * dim + k];
        }

        NSV[1] += -beta * solT_gss * phiV[i];
        NSVold[1] += -beta * solTold_gss * phiV[i];

        for(unsigned  k = 0; k < dim; k++) {
          aResV[k][i] += (- (solV_gss[k] - solVold_gss[k]) * phiV[i] / dt - 0.5 * (NSV[k] + NSVold[k])) * weight;
        }
      } // end phiV_i loop

      // *** phiP_i loop ***
      for(unsigned i = 0; i < nDofsP; i++) {
        for(int k = 0; k < dim; k++) {
          aResP[i] += - (gradSolV_gss[k][k]) * phiP[i]  * weight;
        }
      } // end phiP_i loop

    } // end gauss point loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector

    //copy the value of the adept::adoube aRes in double Res and store them in RES
    Res.resize(nDofsTVP);    //resize

    for(int i = 0; i < nDofsT; i++) {
      Res[i] = -aResT[i].value();
    }

    for(int i = 0; i < nDofsV; i++) {
      for(unsigned  k = 0; k < dim; k++) {
        Res[ i + nDofsT + k * nDofsV ] = -aResV[k][i].value();
      }
    }

    for(int i = 0; i < nDofsP; i++) {
      Res[ i + nDofsT + dim * nDofsV ] = -aResP[i].value();
    }

    RES->add_vector_blocked(Res, sysDof);

    //Extarct and store the Jacobian
    if(assembleMatrix) {
      Jac.resize(nDofsTVP * nDofsTVP);
      // define the dependent variables
      s.dependent(&aResT[0], nDofsT);

      for(unsigned  k = 0; k < dim; k++) {
        s.dependent(&aResV[k][0], nDofsV);
      }

      s.dependent(&aResP[0], nDofsP);

      // define the independent variables
      s.independent(&solT[0], nDofsT);

      for(unsigned  k = 0; k < dim; k++) {
        s.independent(&solV[k][0], nDofsV);
      }

      s.independent(&solP[0], nDofsP);

      // get the and store jacobian matrix (row-major)
      s.jacobian(&Jac[0] , true);
      KK->add_matrix_blocked(Jac, sysDof, sysDof);

      s.clear_independents();
      s.clear_dependents();
    }
  } //end element loop for each process

  RES->close();

  if(assembleMatrix) {
    KK->close();
  }

  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
	// *** get gauss point weight, test function and test function partial derivatives ***
	msh->_finiteElement[kelGeom][soluType]->Jacobian(x,ig,weight,phi,phi_x,phi_xx);
	
	// evaluate the solution, the solution derivatives and the coordinates in the gauss point
	adept::adouble soluGauss = 0;
	vector < adept::adouble > soluGauss_x(dim,0.);
	
	adept::adouble solvGauss = 0;
	vector < adept::adouble > solvGauss_x(dim,0.);

	vector < double > xGauss(dim,0.);
	
	for(unsigned i=0; i<nDofs; i++) {
	  soluGauss+=phi[i]*solu[i];
	  solvGauss+=phi[i]*solv[i];
	  for(unsigned jdim=0; jdim<dim; jdim++) {
	    soluGauss_x[jdim] += phi_x[i*dim+jdim]*solu[i];
	    solvGauss_x[jdim] += phi_x[i*dim+jdim]*solv[i];
	    xGauss[jdim] += x[jdim][i]*phi[i]; 
	  }
	}
        
        double c=.1;
        double Id[2][2]={{1.,0.},{0.,1.}};
	adept::adouble A2 = 1.;
	vector < vector < adept::adouble> > B(dim);
	for(unsigned jdim=0; jdim<dim; jdim++){
	  B[jdim].resize(dim);
	  A2 += soluGauss_x[jdim]*soluGauss_x[jdim]; 
	}
	adept::adouble A = sqrt(A2);
		
	for(unsigned jdim=0; jdim<dim; jdim++){
	  for(unsigned kdim=0; kdim<dim; kdim++){
	     B[jdim][kdim]= Id[jdim][kdim] - (soluGauss_x[jdim] * soluGauss_x[kdim]) / A2;
	  }
	}
		
        // *** phi_i loop ***
	for(unsigned i=0; i<nDofs; i++) {
	 			 
	  adept::adouble nonLinearLaplaceU = 0.;
	  adept::adouble nonLinearLaplaceV = 0.;
	  
	  for(unsigned jdim=0; jdim<dim; jdim++) {
	    
	    nonLinearLaplaceU +=  - 1. / A  * soluGauss_x[jdim] * phi_x[i*dim+jdim]; 
	    
	    nonLinearLaplaceV +=    1. / A * ( - ( B[jdim][0] * solvGauss_x[0] + B[jdim][1] * solvGauss_x[1] ) * phi_x[i*dim+jdim]
					       + ( solvGauss*solvGauss/A2 + c ) * soluGauss_x[jdim] * phi_x[i*dim+jdim] );
	    
	  }
	  	  
	  aResu[i]+= ( 2.*solvGauss/A * phi[i] - nonLinearLaplaceU ) * weight;
	  aResv[i]+= nonLinearLaplaceV * weight;
	  
	} // end phi_i loop
      } // end gauss point loop
    } // endif single element not refined or fine grid loop

    //--------------------------------------------------------------------------------------------------------
    // Add the local Matrix/Vector into the global Matrix/Vector
    
    //copy the value of the adept::adoube aRes in double Res and store
    for(int i=0; i<nDofs; i++) {
      Res[i]       = aResu[i].value();
      Res[nDofs+i] = aResv[i].value();
    }
    RES->add_vector_blocked(Res,KKDof);
    
    if(assembleMatrix) {
      
      // define the dependent variables
      s.dependent(&aResu[0], nDofs);
      s.dependent(&aResv[0], nDofs);
      
      // define the independent variables
      s.independent(&solu[0], nDofs);   
      s.independent(&solv[0], nDofs);   
      // get the jacobian matrix (ordered by column)
      s.jacobian(&Jac[0]);	
     
      // get the jacobian matrix (ordered by raw, i.e. Jact=Jac^t)
      for (int inode=0;inode<2*nDofs;inode++){
	for (int jnode=0;jnode<2*nDofs;jnode++){
	   Jact[inode*2*nDofs+jnode]=-Jac[jnode*2*nDofs+inode];
	}
      } 
      //store Jact in the global matrix KK
      KK->add_matrix_blocked(Jact,KKDof,KKDof);
      
      s.clear_independents();
      s.clear_dependents();
    }
  } //end element loop for each process

  RES->close();
  if( assembleMatrix ) KK->close();
  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
    
      for(int k = 0; k < n_unknowns; k++) {
      L2G_dofmap[k].resize(Sol_n_el_dofs[k]);

      Res_el[SolPdeIndex[k]].resize(Sol_n_el_dofs[k]);
      memset(& Res_el[SolPdeIndex[k]][0], 0., Sol_n_el_dofs[k] * sizeof(double) );


      Jac_el[SolPdeIndex[k]][SolPdeIndex[k]].resize(Sol_n_el_dofs[k] * Sol_n_el_dofs[k]);
      memset(& Jac_el[SolPdeIndex[k]][SolPdeIndex[k]][0], 0., Sol_n_el_dofs[k] * Sol_n_el_dofs[k] * sizeof(double));
    }

  //all vars###################################################################    

      // *** quadrature loop ***
      for(unsigned ig = 0; ig < ml_prob.GetQuadratureRule(ielGeom).GetGaussPointsNumber(); ig++) {
          
      // *** get gauss point weight, test function and test function partial derivatives ***
      for(int fe=0; fe < NFE_FAMS; fe++) {
         msh->_finiteElement[ielGeom][fe]->Jacobian(coords_at_dofs,ig,weight_qp,phi_fe_qp[fe],phi_x_fe_qp[fe],phi_xx_fe_qp[fe]);
      }
      //HAVE TO RECALL IT TO HAVE BIQUADRATIC JACOBIAN
         msh->_finiteElement[ielGeom][coords_fe_type]->Jacobian(coords_at_dofs,ig,weight_qp,phi_fe_qp[coords_fe_type],phi_x_fe_qp[coords_fe_type],phi_xx_fe_qp[coords_fe_type]);

   //========= fill gauss value quantities ==================   
   std::fill(sol_qp.begin(), sol_qp.end(), 0.);
   std::fill(sol_old_qp.begin(), sol_old_qp.end(), 0.);
   for (unsigned  k = 0; k < n_unknowns; k++) { std::fill(sol_grad_qp[k].begin(), sol_grad_qp[k].end(), 0.); 
                                                std::fill(sol_old_grad_qp[k].begin(), sol_old_grad_qp[k].end(), 0.);
                                            }
    
    for (unsigned  k = 0; k < n_unknowns; k++) {
	for (unsigned i = 0; i < Sol_n_el_dofs[k]; i++) {
	                                                         sol_qp[k]    +=     sol_eldofs[k][i] *   phi_fe_qp[SolFEType[k]][i];
	                                                     sol_old_qp[k]    += sol_old_eldofs[k][i] *   phi_fe_qp[SolFEType[k]][i];
                   for (unsigned d = 0; d < dim; d++) {      sol_grad_qp[k][d] +=     sol_eldofs[k][i] * phi_x_fe_qp[SolFEType[k]][i * dim + d]; 
                                                         sol_old_grad_qp[k][d] += sol_old_eldofs[k][i] * phi_x_fe_qp[SolFEType[k]][i * dim + d]; 
                                     }
       }        
    }
 
  //========= fill gauss value quantities ==================
         
	// *** phi_i loop ***
	for(unsigned i = 0; i < nDof_max; i++) {

    //BEGIN RESIDUALS A block ===========================
	    double Lap_rhs_i = 0.;
	    double Lap_old_rhs_i = 0.;
	    for(unsigned d = 0;  d < dim;  d++) {
	      Lap_rhs_i     += phi_x_fe_qp[SolFEType[0]][i * dim + d] *     sol_grad_qp[0][d];
	      Lap_old_rhs_i += phi_x_fe_qp[SolFEType[0]][i * dim + d] * sol_old_grad_qp[0][d];
	    }
	    Res_el[SolPdeIndex[0]][i] +=  - weight_qp * ( 
                       delta_g_term * Singularity::function(sol_qp[0]) * Singularity::g_vc(sol_qp[0])     * phi_fe_qp[ SolFEType[0] ][i]
                     - delta_g_term * Singularity::function(sol_qp[0]) * Singularity::g_vc(sol_old_qp[0]) * phi_fe_qp[ SolFEType[0] ][i]
                     + deltat_term * Singularity::function(sol_qp[0]) * dt                               * phi_fe_qp[ SolFEType[0] ][i]
                     + lapl_term *  dt * Lap_rhs_i
            );
    //END RESIDUALS A block ===========================


	    // *** phi_j loop ***
	    for(unsigned j = 0; j<nDof_max; j++) {
            
          double Lap_mat_i_j = 0.;
          for(unsigned d = 0; d < dim; d++) Lap_mat_i_j += phi_x_fe_qp[SolFEType[0]][i * dim + d] *
                                                           phi_x_fe_qp[SolFEType[0]][j * dim + d];


          Jac_el[SolPdeIndex[0]][SolPdeIndex[0]][i * Sol_n_el_dofs[0] + j] += weight_qp * (
                    + delta_g_term * phi_fe_qp[ SolFEType[0] ][i] * Singularity::derivative(  phi_fe_qp[ SolFEType[0] ][j] ) * phi_fe_qp[ SolFEType[0] ][j] * Singularity::g_vc( phi_fe_qp[ SolFEType[0] ][j] ) 
                    + delta_g_term * phi_fe_qp[ SolFEType[0] ][i] * Singularity::function  (  phi_fe_qp[ SolFEType[0] ][j] ) * Singularity::g_vc_derivative( phi_fe_qp[ SolFEType[0] ][j] ) * phi_fe_qp[ SolFEType[0] ][j] 
                    
                    - delta_g_term * phi_fe_qp[ SolFEType[0] ][i] * Singularity::derivative(  phi_fe_qp[ SolFEType[0] ][j] ) * phi_fe_qp[ SolFEType[0] ][j] * Singularity::g_vc( sol_old_qp[0] ) 
                    
                    + deltat_term * phi_fe_qp[ SolFEType[0] ][i] * Singularity::derivative(  phi_fe_qp[ SolFEType[0] ][j] ) * phi_fe_qp[ SolFEType[0] ][j] * dt 
                    + lapl_term * dt * Lap_mat_i_j  
                                                                                          );
	      
  	    }    //end phij loop

        
    }        //end phii loop
    
  }          // end gauss point loop

//--------------------------------------------------------------------------------------------------------
    //Sum the local matrices/vectors into the Global Matrix/Vector
    for(unsigned ivar=0; ivar < n_unknowns; ivar++) {
      RES->add_vector_blocked(Res_el[SolPdeIndex[ivar]],L2G_dofmap[ivar]);
      JAC->add_matrix_blocked(Jac_el[SolPdeIndex[ivar]][SolPdeIndex[ivar]],L2G_dofmap[ivar],L2G_dofmap[ivar]);
    }
    //--------------------------------------------------------------------------------------------------------
  } //end list of elements loop for each subdomain

  JAC->close();
  RES->close();
  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........

	      for(unsigned ivar=0; ivar<dim; ivar++) {    
		B[SolPdeIndex[ivar]][SolPdeIndex[ivar]][i*nve2+j] += IRe*Lap+ NavierStokes*newton*Adv1;
		if(nwtn_alg==2){
		  // Advection term II
		  B[SolPdeIndex[ivar]][SolPdeIndex[ivar]][i*nve2+j]       += Adv2*gradSolVAR[ivar][ivar];
		  for(unsigned ivar2=1; ivar2<dim; ivar2++) {
		    B[SolPdeIndex[ivar]][SolPdeIndex[(ivar+ivar2)%dim]][i*nve2+j] += Adv2*gradSolVAR[ivar][(ivar+ivar2)%dim];
		  }
		}
	      }
  	    } //end phij loop
	    
	    // *** phi1_j loop ***
	    for(unsigned j=0; j<nve1; j++){
	      for(unsigned ivar=0; ivar<dim; ivar++) {
		B[SolPdeIndex[ivar]][SolPdeIndex[dim]][i*nve1+j] -= gradphi2[i*dim+ivar]*phi1[j]*Weight2;
	      }
	    } //end phi1_j loop
	  } // endif assemble_matrix
	} //end phii loop
  

	// *** phi1_i loop ***
	for(unsigned i=0; i<nve1; i++){
	  //BEGIN RESIDUALS B block ===========================
	  double div = 0;
	  for(unsigned ivar=0; ivar<dim; ivar++) {
	    div += gradSolVAR[ivar][ivar];
	  }
	  F[SolPdeIndex[dim]][i]+= (phi1[i]*div +penalty*ILambda*phi1[i]*SolVAR[dim])*Weight2;
	  //END RESIDUALS  B block ===========================
	  
	  if(assemble_matrix){
	    // *** phi_j loop ***
	    for(unsigned j=0; j<nve2; j++) {
	      for(unsigned ivar=0; ivar<dim; ivar++) {
		B[SolPdeIndex[dim]][SolPdeIndex[ivar]][i*nve2+j]-= phi1[i]*gradphi2[j*dim+ivar]*Weight2;
	      }
	    }  //end phij loop
	  } // endif assemble_matrix
	}  //end phi1_i loop
	
	if(assemble_matrix * penalty){  //block nve1 nve1
	  // *** phi_i loop ***
	  for(unsigned i=0; i<nve1; i++){
	    // *** phi_j loop ***
	    for(unsigned j=0; j<nve1; j++){
	      B[SolPdeIndex[dim]][SolPdeIndex[dim]][i*nve1+j]-= ILambda*phi1[i]*phi1[j]*Weight2;
	    }
	  }
	}   //end if penalty
      }  // end gauss point loop
      
      //--------------------------------------------------------------------------------------------------------
      // Boundary Integral --> to be addded
      //number of faces for each type of element
//       if (igrid==gridn || !myel->GetRefinedElementIndex(kel) ) {
//      
// 	unsigned nfaces = myel->GetElementFaceNumber(kel);
// 
// 	// loop on faces
// 	for(unsigned jface=0;jface<nfaces;jface++){ 
// 	  
// 	  // look for boundary faces
// 	  if(myel->GetFaceElementIndex(kel,jface)<0){
// 	    for(unsigned ivar=0; ivar<dim; ivar++) {
// 	      ml_prob.ComputeBdIntegral(pdename, &Solname[ivar][0], kel, jface, level, ivar);
// 	    }
// 	  }
// 	}	
//       }
      //--------------------------------------------------------------------------------------------------------
    } // endif single element not refined or fine grid loop
    //--------------------------------------------------------------------------------------------------------
    //Sum the local matrices/vectors into the Global Matrix/Vector
    for(unsigned ivar=0; ivar<dim; ivar++) {
      myRES->add_vector_blocked(F[SolPdeIndex[ivar]],KK_dof[ivar]);
      if(assemble_matrix){
	myKK->add_matrix_blocked(B[SolPdeIndex[ivar]][SolPdeIndex[ivar]],KK_dof[ivar],KK_dof[ivar]);  
	myKK->add_matrix_blocked(B[SolPdeIndex[ivar]][SolPdeIndex[dim]],KK_dof[ivar],KK_dof[dim]);
	myKK->add_matrix_blocked(B[SolPdeIndex[dim]][SolPdeIndex[ivar]],KK_dof[dim],KK_dof[ivar]);
	if(nwtn_alg==2){
	  for(unsigned ivar2=1; ivar2<dim; ivar2++) {
	    myKK->add_matrix_blocked(B[SolPdeIndex[ivar]][SolPdeIndex[(ivar+ivar2)%dim]],KK_dof[ivar],KK_dof[(ivar+ivar2)%dim]);  
	  }
	}
      }
    }
    //Penalty
    if(assemble_matrix*penalty) myKK->add_matrix_blocked(B[SolPdeIndex[dim]][SolPdeIndex[dim]],KK_dof[dim],KK_dof[dim]);
    myRES->add_vector_blocked(F[SolPdeIndex[dim]],KK_dof[dim]);
    //--------------------------------------------------------------------------------------------------------  
  } //end list of elements loop for each subdomain
  
  
  if(assemble_matrix) myKK->close();
  myRES->close();
  // ***************** END ASSEMBLY *******************
}

//.........这里部分代码省略.........
	      //=========== delta_adjoint row ===========================
		Jac[ assemble_jacobian<double,double>::jac_row_col_index(Sol_n_el_dofs, nDof_AllVars, pos_adj, pos_state, i, j) ]  += m_b_f[pos_adj][pos_state] * weight_qp * (
                                                                                                                                      assemble_jacobian<double,double>::laplacian_row_col(phi_x_fe_qp_vec,phi_x_fe_qp_vec, pos_adj, pos_state, i, j, dim)
                                                                                                                                      - phi_fe_qp[pos_adj][i] *
                                                                                                                                        nonlin_term_derivative(phi_fe_qp[pos_state][j]) * 
                                                                                                                                        sol_qp[pos_ctrl]
                                                                                                                                    );

        Jac[ assemble_jacobian<double,double>::jac_row_col_index(Sol_n_el_dofs, nDof_AllVars, pos_adj, pos_ctrl, i, j)  ]  += m_b_f[pos_adj][pos_ctrl] * weight_qp * ( 
                                                                                                                                      - phi_fe_qp[pos_adj][i] * 
                                                                                                                                        phi_fe_qp[pos_ctrl][j] * 
                                                                                                                                        nonlin_term_function(sol_qp[pos_state])
                                                                                                                                   ); 
		     
	      
            } // end phi_j loop
          } // endif assemble_matrix

        } // end phi_i loop
        
      } // end gauss point loop

      
      
    //std::vector<double> Res_ctrl (Sol_n_el_dofs[pos_ctrl]); std::fill(Res_ctrl.begin(),Res_ctrl.end(), 0.);
    for (unsigned i = 0; i < sol_eldofs[pos_ctrl].size(); i++) {
       unsigned n_els_that_node = 1;
     if ( control_el_flag == 1) {
// 	Res[Sol_n_el_dofs[pos_state] + i] += - ( + n_els_that_node * ineq_flag * sol_eldofs[pos_mu][i] /*- ( 0.4 + sin(M_PI * x[0][i]) * sin(M_PI * x[1][i]) )*/ );
// 	Res_ctrl[i] =  Res[Sol_n_el_dofs[pos_state] + i];
      }
    }
    

 //========== end of integral-based part

                          RES->add_vector_blocked(Res, L2G_dofmap_AllVars);
      if (assembleMatrix)  KK->add_matrix_blocked(Jac, L2G_dofmap_AllVars, L2G_dofmap_AllVars);
    
    
 //========== dof-based part, without summation
 
 //============= delta_mu row ===============================
      std::vector<double> Res_mu (Sol_n_el_dofs[pos_mu]); std::fill(Res_mu.begin(),Res_mu.end(), 0.);
      
    for (unsigned i = 0; i < sol_actflag.size(); i++) {
      if (sol_actflag[i] == 0) {  //inactive
         Res_mu [i] = (- ineq_flag) * ( 1. * sol_eldofs[pos_mu][i] ); 
// 	 Res_mu [i] = Res[res_row_index(Sol_n_el_dofs,pos_mu,i)]; 
      }
      else if (sol_actflag[i] == 1) {  //active_a 
	     Res_mu [i] = (- ineq_flag) * c_compl * ( sol_eldofs[pos_ctrl][i] - ctrl_lower[i]);
      }
      else if (sol_actflag[i] == 2) {  //active_b 
	     Res_mu [i] = (- ineq_flag) * c_compl * ( sol_eldofs[pos_ctrl][i] - ctrl_upper[i]);
      }
    }
    
    RES->insert(Res_mu, L2G_dofmap[pos_mu]);
    
 //============= delta_ctrl - mu ===============================
 KK->matrix_set_off_diagonal_values_blocked(L2G_dofmap[pos_ctrl], L2G_dofmap[pos_mu], m_b_f[pos_ctrl][pos_mu] * ineq_flag * 1.);
  
 //============= delta_mu - ctrl row ===============================
 for (unsigned i = 0; i < sol_actflag.size(); i++) if (sol_actflag[i] != 0 ) sol_actflag[i] = m_b_f[pos_mu][pos_ctrl] * ineq_flag * c_compl;    
  
 KK->matrix_set_off_diagonal_values_blocked(L2G_dofmap[pos_mu], L2G_dofmap[pos_ctrl], sol_actflag);

 //============= delta_mu - mu row ===============================
  for (unsigned i = 0; i < sol_actflag.size(); i++) sol_actflag[i] =  m_b_f[pos_mu][pos_mu] * ( ineq_flag * (1 - sol_actflag[i]/c_compl)  + (1-ineq_flag) * 1. ); //can do better to avoid division, maybe use modulo operator 

  KK->matrix_set_off_diagonal_values_blocked(L2G_dofmap[pos_mu], L2G_dofmap[pos_mu], sol_actflag );
  
  } //end element loop for each process
  
  RES->close();

  if (assembleMatrix) KK->close();
//  std::ostringstream mat_out; mat_out << ml_prob.GetFilesHandler()->GetOutputPath() << "/" << "jacobian" << mlPdeSys->GetNonlinearIt()  << ".txt";
//   KK->print_matlab(mat_out.str(),"ascii"); 
//    KK->print();
  
  // ***************** END ASSEMBLY *******************
  
  unsigned int global_ctrl_size = pdeSys->KKoffset[pos_ctrl+1][iproc] - pdeSys->KKoffset[pos_ctrl][iproc];
  
  std::vector<double>  one_times_mu(global_ctrl_size, 0.);
  std::vector<int>        positions(global_ctrl_size);
//  double position_mu_i;
  for (unsigned i = 0; i < positions.size(); i++) {
    positions[i] = pdeSys->KKoffset[pos_ctrl][iproc] + i;
//     position_mu_i = pdeSys->KKoffset[pos_mu][iproc] + i;
//     std::cout << position_mu_i << std::endl;
    one_times_mu[i] =  m_b_f[pos_ctrl][pos_mu] * ineq_flag * 1. * (*sol->_Sol[SolIndex[pos_mu]])(i/*position_mu_i*/) ;
  }
    RES->add_vector_blocked(one_times_mu, positions);
//     RES->print();
    
  return;
}