C++ VectorXi::rows方法代码示例

void drawConstraints(igl::viewer::Viewer &viewer)
{
  for (int n=0; n<2; ++n)
  {
    Eigen::MatrixXd Bc = igl::slice(B, b, 1);
    Eigen::MatrixXd color;
    color.setZero(b.rows(),3);
    color.col(2).setOnes();
    for (int i =0; i<b.rows(); ++i)
      if (blevel[i] ==1 && n>0)
        color.row(i)<<0.7,0.7,0.7;
    // Eigen::RowVector3d color; color<<0.5,0.5,0.5;
    viewer.data.add_edges(Bc - global_scale*bc.block(0,n*3,bc.rows(),3), Bc + global_scale*bc.block(0,n*3,bc.rows(),3) , color);
  }

}

bool update_arap()
{
  using namespace Eigen;
  using namespace igl;
  using namespace std;
  MatrixXd bc(num_in_selection(S),V.cols());
  // get b from S
  {
    if(!paused)
    {
      t = get_seconds()-pause_time;
    }
    int bi = 0;
    for(int v = 0;v<S.rows(); v++)
    {
      if(S(v) >= 0)
      {
        bc.row(bi) = V.row(v);
        switch(S(v))
        {
          case 0:
          {
            const double r = mid(0)*0.25;
            bc(bi,0) += r*cos(0.5*t*2.*PI);
            bc(bi,1) -= r+r*sin(0.5*t*2.*PI);
            break;
          }
          case 1:
          {
            const double r = mid(1)*0.15;
            bc(bi,1) += r+r*cos(0.15*t*2.*PI);
            bc(bi,2) -= r*sin(0.15*t*2.*PI);

            //// Pull-up
            //bc(bi,0) += 0.42;//mid(0)*0.5;
            //bc(bi,1) += 0.55;//mid(0)*0.5;
            //// Bend
            //Vector3d t(-1,0,0);
            //Quaterniond q(AngleAxisd(PI/1.5,Vector3d(0,1.0,0.1).normalized()));
            //const Vector3d a = bc.row(bi);
            //bc.row(bi) = (q*(a-t) + t) + Vector3d(1.5,0.1,0.9);


            break;
          }
          default:
          break;
        }
        bi++;
      }
    }
  }
  if(!arap_solve(bc,arap_data,U))
  {
    cerr<<"arap_solve failed."<<endl;
    return false;
  }
  per_face_normals(U,F,N);
  return true;
}

void NuTo::Structure::ElementCreate(int rElementNumber, int rInterpolationTypeId, const Eigen::VectorXi& rNodeNumbers,
                                    const Eigen::MatrixXd& rKnots, const Eigen::VectorXi& rKnotIDs)
{
    // convert node numbers to pointer
    std::vector<NodeBase*> nodeVector;
    for (int iNode = 0; iNode < rNodeNumbers.rows(); iNode++)
        nodeVector.push_back(NodeGetNodePtr(rNodeNumbers(iNode)));

    boost::ptr_map<int, InterpolationType>::iterator itIterator = mInterpolationTypeMap.find(rInterpolationTypeId);

    if (itIterator == mInterpolationTypeMap.end())
        throw NuTo::Exception(__PRETTY_FUNCTION__, "Interpolation type does not exist.");

    InterpolationType& interpolationType = *itIterator->second;

    if (not interpolationType.IsDof(Node::eDof::COORDINATES))
        throw NuTo::Exception(__PRETTY_FUNCTION__, "COORDINATE interpolation required.");

    unsigned int numNodesCoordinates = interpolationType.Get(Node::eDof::COORDINATES).GetNumNodes();
    if (numNodesCoordinates != nodeVector.size())
        throw NuTo::Exception(__PRETTY_FUNCTION__, "COORDINATE interpolation requires " +
                                                           std::to_string(numNodesCoordinates) + " nodes. " +
                                                           std::to_string(nodeVector.size()) + " are provided.");

    interpolationType.ClearCache();
    const auto& integrationType = *GetPtrIntegrationType(interpolationType.GetStandardIntegrationType());

    ElementBase* ptrElement = nullptr;
    switch (interpolationType.GetShapeType())
    {
    case NuTo::Interpolation::eShapeType::SPRING:
        throw NuTo::Exception(__PRETTY_FUNCTION__, "Element1DSpring currently not implemented.");
        break;
    case NuTo::Interpolation::eShapeType::TRUSS1D:
    case NuTo::Interpolation::eShapeType::TRUSSXD:
    case NuTo::Interpolation::eShapeType::TRIANGLE2D:
    case NuTo::Interpolation::eShapeType::QUAD2D:
    case NuTo::Interpolation::eShapeType::TETRAHEDRON3D:
    case NuTo::Interpolation::eShapeType::BRICK3D:
    case NuTo::Interpolation::eShapeType::INTERFACE:
        throw NuTo::Exception(__PRETTY_FUNCTION__,
                              "Please use approriate functions for element creation, this is IGA implementation.");
        break;
    case NuTo::Interpolation::eShapeType::IGA1D:
        ptrElement = new ContinuumElementIGA<1>(nodeVector, rKnots, rKnotIDs, interpolationType, integrationType,
                                                GetDofStatus());
        ptrElement->CheckElement();
        break;
    case NuTo::Interpolation::eShapeType::IGA2D:
        ptrElement = new ContinuumElementIGA<2>(nodeVector, rKnots, rKnotIDs, interpolationType, integrationType,
                                                GetDofStatus());
        ptrElement->CheckElement();
        break;
    default:
        throw Exception(__PRETTY_FUNCTION__, "invalid dimension.");
    }

    mElementMap.insert(rElementNumber, ptrElement);
}

int num_in_selection(const Eigen::VectorXi & S)
{
  int count = 0;
  for(int v = 0;v<S.rows(); v++)
  {
    if(S(v) >= 0)
    {
      count++;
    }
  }
  return count;
}

void randomly_color(
  const Eigen::VectorXi & CC,
  Eigen::MatrixXd & C)
{
  using namespace Eigen;
  using namespace igl;
  using namespace std;
  VectorXi I;
  srand ( unsigned ( time(0) ) );
  double num_cc = (double)CC.maxCoeff()+1.0;
  randperm(num_cc,I);
  C.resize(CC.rows(),3);
  for(int f = 0;f<CC.rows();f++)
  {
    jet(
      (double)I(CC(f))/num_cc,
      C(f,0),
      C(f,1),
      C(f,2));
  }
}

		int findValue(Eigen::VectorXi vec, int value) {
			int toReturn  = -1;
			int currNdx = 0;
			while (toReturn == -1 && currNdx < vec.rows()) {
				if (vec(currNdx) == value) {
					toReturn = currNdx;
				} else {
					currNdx++;
				}
			}
			return toReturn;
		}

        void addEdgeType( CEdgeTypePtr edgeType, Eigen::VectorXi &typeOfEdgeFeatures )
        {
            // Consistency checks
            assert( edgeType->getNumberOfFeatures() == typeOfEdgeFeatures.rows() );

            ITERATE_SIZE_T(typeOfEdgeFeatures)
                assert( typeOfEdgeFeatures(i) <= 3 );

            // Add the edge type to the vector of edge types
            m_edgeTypes.push_back( edgeType );
            m_typesOfEdgeFeatures[ edgeType ] = typeOfEdgeFeatures;
        }

void parse_rhs(
  const int nrhs, 
  const mxArray *prhs[], 
  Eigen::MatrixXd & V,
  Eigen::MatrixXi & Ele,
  Eigen::MatrixXd & Q,
  Eigen::MatrixXd & bb_mins,
  Eigen::MatrixXd & bb_maxs,
  Eigen::VectorXi & elements)
{
  using namespace std;
  using namespace igl;
  using namespace igl::matlab;
  mexErrMsgTxt(nrhs >= 3, "The number of input arguments must be >=3.");

  const int dim = mxGetN(prhs[0]);
  mexErrMsgTxt(dim == 3 || dim == 2,
    "Mesh vertex list must be #V by 2 or 3 list of vertex positions");

  mexErrMsgTxt(dim+1 == mxGetN(prhs[1]),
    "Mesh \"face\" simplex size must equal dimension+1");

  parse_rhs_double(prhs,V);
  parse_rhs_index(prhs+1,Ele);
  parse_rhs_double(prhs+2,Q);
  mexErrMsgTxt(Q.cols() == dim,"Dimension of Q should match V");
  if(nrhs > 3)
  {
    mexErrMsgTxt(nrhs >= 6, "The number of input arguments must be 3 or >=6.");
    parse_rhs_double(prhs+3,bb_mins);
    if(bb_mins.size()>0)
    {
      mexErrMsgTxt(bb_mins.cols() == dim,"Dimension of bb_mins should match V");
      mexErrMsgTxt(bb_mins.rows() >= Ele.rows(),"|bb_mins| should be > |Ele|");
    }
    parse_rhs_double(prhs+4,bb_maxs);
    mexErrMsgTxt(bb_maxs.cols() == bb_mins.cols(),
      "|bb_maxs| should match |bb_mins|");
    mexErrMsgTxt(bb_mins.rows() == bb_maxs.rows(),
      "|bb_mins| should match |bb_maxs|");
    parse_rhs_index(prhs+5,elements);
    mexErrMsgTxt(elements.cols() == 1,"Elements should be column vector");
    mexErrMsgTxt(bb_mins.rows() == elements.rows(),
      "|bb_mins| should match |elements|");
  }else
  {
    // Defaults
    bb_mins.resize(0,dim);
    bb_maxs.resize(0,dim);
    elements.resize(0,1);
  }
}

IGL_INLINE void igl::forward_kinematics(
  const Eigen::MatrixXd & C,
  const Eigen::MatrixXi & BE,
  const Eigen::VectorXi & P,
  const std::vector<
    Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
  std::vector<
    Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > & vQ,
  std::vector<Eigen::Vector3d> & vT)
{
  using namespace std;
  using namespace Eigen;
  const int m = BE.rows(); 
  assert(m == P.rows());
  assert(m == (int)dQ.size());
  vector<bool> computed(m,false);
  vQ.resize(m);
  vT.resize(m);
  function<void (int) > fk_helper = [&] (int b)
  {
    if(!computed[b])
    {
      if(P(b) < 0)
      {
        // base case for roots
        vQ[b] = dQ[b];
        const Vector3d r = C.row(BE(b,0)).transpose();
        vT[b] = r-dQ[b]*r;
      }else
      {
        // Otherwise first compute parent's
        const int p = P(b);
        fk_helper(p);
        vQ[b] = vQ[p] * dQ[b];
        const Vector3d r = C.row(BE(b,0)).transpose();
        vT[b] = vT[p] - vQ[b]*r + vQ[p]*r;
      }
      computed[b] = true;
    }
  };
  for(int b = 0;b<m;b++)
  {
    fk_helper(b);
  }
}

IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
                             const Eigen::MatrixXi &F,
                             const Eigen::VectorXi& b,
                             const Eigen::MatrixXd& bc,
                             const Eigen::VectorXi &I,
                             Eigen::MatrixXd &output)
{
  Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
  Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);

  for(unsigned i=0; i<b.size();++i)
  {
    isConstrained(b(i)) = 1;
    cfW.row(b(i)) << bc.row(i);
  }
  int n = I.rows();
  igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
  pvff.solve(isConstrained, cfW, I, output);
}

bool
writeLabelToFile (const std::string &file, const Eigen::VectorXi &vector)
{
    std::ofstream out (file.c_str ());
    if (!out)
    {
        std::cout << "Cannot open file.\n";
        return false;
    }
    size_t i ,j;

    for (i = 0; i < vector.rows(); i++)
    {

            out << vector (i);


        out<<" ";
    }
    out.close ();

    return true;
}

int main(int argc, char *argv[])
{
  using namespace std;
  // Load a mesh in OFF format
  igl::readOBJ(TUTORIAL_SHARED_PATH "/camel_b.obj", V, F);

  Eigen::MatrixXd bnd_uv, uv_init;

  Eigen::VectorXd M;
  igl::doublearea(V, F, M);
  std::vector<std::vector<int>> all_bnds;
  igl::boundary_loop(F, all_bnds);

  // Heuristic primary boundary choice: longest
  auto primary_bnd = std::max_element(all_bnds.begin(), all_bnds.end(), [](const std::vector<int> &a, const std::vector<int> &b) { return a.size()<b.size(); });

  Eigen::VectorXi bnd = Eigen::Map<Eigen::VectorXi>(primary_bnd->data(), primary_bnd->size());

  igl::map_vertices_to_circle(V, bnd, bnd_uv);
  bnd_uv *= sqrt(M.sum() / (2 * igl::PI));
  if (all_bnds.size() == 1)
  {
    if (bnd.rows() == V.rows()) // case: all vertex on boundary
    {
      uv_init.resize(V.rows(), 2);
      for (int i = 0; i < bnd.rows(); i++)
        uv_init.row(bnd(i)) = bnd_uv.row(i);
    }
    else
    {
      igl::harmonic(V, F, bnd, bnd_uv, 1, uv_init);
      if (igl::flipped_triangles(uv_init, F).size() != 0)
        igl::harmonic(F, bnd, bnd_uv, 1, uv_init); // fallback uniform laplacian
    }
  }
  else
  {
    // if there is a hole, fill it and erase additional vertices.
    all_bnds.erase(primary_bnd);
    Eigen::MatrixXi F_filled;
    igl::topological_hole_fill(F, bnd, all_bnds, F_filled);
    igl::harmonic(F_filled, bnd, bnd_uv ,1, uv_init);
    uv_init = uv_init.topRows(V.rows());
  }

  Eigen::VectorXi b; Eigen::MatrixXd bc;
  igl::scaf_precompute(V, F, uv_init, scaf_data, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0);

  // Plot the mesh
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V, F);
  const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows());
  viewer.data().set_uv(V_uv);
  viewer.callback_key_down = &key_down;

  // Enable wireframe
  viewer.data().show_lines = true;

  // Draw checkerboard texture
  viewer.data().show_texture = true;


  std::cerr << "Press space for running an iteration." << std::endl;
  std::cerr << "Press 1 for Mesh 2 for UV" << std::endl;

  // Launch the viewer
  viewer.launch();
}

void SliceStack::computeLaplace(const Eigen::MatrixXd &TV, const Eigen::MatrixXi &TT,
                                const Eigen::MatrixXi &TF, const Eigen::VectorXi &TO,
                                Eigen::VectorXd &Z, const set<int> &allowed) {
  bool laplace_DEBUG = false;
	Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
	Eigen::IOFormat LongFmt(10, 0, ", ", "\n", "[", "]");
	Eigen::IOFormat RFmt(4, 0, ", ", ", ", "", "", "(", ")");

	assert(TO.rows() == TV.rows());

	std::vector<int> known_v;
	std::vector<double> known_c_v;
  auto mx = TV.colwise().maxCoeff();
  auto mn = TV.colwise().minCoeff();

	for (int i = 0; i < TO.rows(); ++i) {
		if (allowed.size() == 0 && TO(i) != GLOBAL::nonoriginal_marker) {
			known_v.push_back(i);
			known_c_v.push_back(GLOBAL::inside_temp);
		}
    else if (allowed.find(TO(i)) != allowed.end()) {
			known_v.push_back(i);
			known_c_v.push_back(GLOBAL::inside_temp);
    }
    else if (TV(i,2) == mx(2) || TV(i,2) == mn(2)) {
			known_v.push_back(i);
			known_c_v.push_back(GLOBAL::outside_temp);
		}
	}

	if (laplace_DEBUG) 
    printf("done! Number of known values is %lu/%lu\n",
           known_v.size(), TV.rows());

	Eigen::VectorXi known(known_v.size());
	Eigen::VectorXd known_c(known_v.size());

	for (int i = 0; i < known_c.size(); ++i) {
		known(i) = known_v[i];
		known_c(i) = known_c_v[i];
	}

	Eigen::SparseMatrix<double> L(TV.rows(), TV.rows());
	// Set non-diag elements to 1 if connected, 0 otherwise
	// Use the tets instead of the faces
	for (int i = 0; i < TT.rows(); ++i) {
		L.coeffRef(TT(i,0), TT(i,1)) = -1; L.coeffRef(TT(i,1), TT(i,0)) = -1;
		L.coeffRef(TT(i,1), TT(i,2)) = -1; L.coeffRef(TT(i,2), TT(i,1)) = -1;
		L.coeffRef(TT(i,2), TT(i,3)) = -1; L.coeffRef(TT(i,3), TT(i,2)) = -1;
		L.coeffRef(TT(i,3), TT(i,0)) = -1; L.coeffRef(TT(i,0), TT(i,3)) = -1;
	}

	// Set diag elements to valence of entry
	for (int i = 0; i < TV.rows(); ++i) {
		L.coeffRef(i,i) = -L.row(i).sum();
	}

  if (laplace_DEBUG) {
    printf("done! Number non-zeros is %ld\n", L.nonZeros());
    printf("Solving energy constraints...");
  }

	// Solve energy constraints.
	igl::min_quad_with_fixed_data<double> mqwf;
	// Linear term is 0
	Eigen::VectorXd B = Eigen::VectorXd::Zero(TV.rows(), 1);
	// Empty Constraints
	Eigen::VectorXd Beq;
	Eigen::SparseMatrix<double> Aeq;

	if (!igl::min_quad_with_fixed_precompute(L, known, Aeq, false, mqwf))
		fprintf(stderr,"ERROR: fixed_precompute didn't work!\n");

	igl::min_quad_with_fixed_solve(mqwf,B,known_c,Beq, Z);

  if (laplace_DEBUG)
    printf("fixed_solve complete.\n");
}

bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
{
  using namespace std;
  using namespace Eigen;

  if (key <'1' || key >'9')
    return false;

  viewer.data.lines.resize(0,9);

  int num = key  - '0';

  // Interpolate
  std::cerr << "Interpolating " << num  << "-PolyVector field" << std::endl;

  VectorXi b(3);
  b << 1511, 603, 506;

  int numConstraintsToGenerate;
  // if it's not a 2-PV or a 1-PV, include a line direction (2 opposite vectors)
  // in the field
  if (num>=5)
    numConstraintsToGenerate  = num-2;
  else
    if (num>=3)
      numConstraintsToGenerate  = num-1;
    else
      numConstraintsToGenerate  = num;


  MatrixXd bc(b.size(),numConstraintsToGenerate*3);
  for (unsigned i=0; i<b.size(); ++i)
  {
    VectorXd t = random_constraints(B1.row(b(i)),B2.row(b(i)),numConstraintsToGenerate);
    bc.row(i) = t;
  }
  VectorXi rootsIndex(num);
  for (int i =0; i<numConstraintsToGenerate; ++i)
    rootsIndex[i] = i+1;
  if (num>=5)
    rootsIndex[num-2] = -2;
    if (num>=3)
      rootsIndex[num-1] = -1;

  // Interpolated PolyVector field
  Eigen::MatrixXd pvf;
  igl::n_polyvector_general(V, F, b, bc, rootsIndex, pvf);

  ofstream ofs;
  ofs.open("pvf.txt", ofstream::out);
  ofs<<pvf;
  ofs.close();
  igl::writeOFF("pvf.off",V,F);
  
  // Highlight in red the constrained faces
  MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  for (unsigned i=0; i<b.size();++i)
    C.row(b(i)) << 1, 0, 0;
  viewer.data.set_colors(C);

  for (int n=0; n<num; ++n)
  {
//    const MatrixXd &VF = pvf.block(0,n*3,F.rows(),3);
    MatrixXd VF = MatrixXd::Zero(F.rows(),3);
    for (unsigned i=0; i<b.size(); ++i)
      VF.row(b[i]) = pvf.block(b[i],n*3,1,3);
    
    for (int i=0; i<samples.rows(); ++i)
      VF.row(samples[i]) = pvf.block(samples[i],n*3,1,3);
    
    VectorXd c = VF.rowwise().norm();
    MatrixXd C2;
    igl::jet(c,1,1+rand_factor,C2);
    // viewer.data.add_edges(B - global_scale*VF, B + global_scale*VF , C2);
    viewer.data.add_edges(B, B + global_scale*VF , C2);
  }


  return false;
}

int NuTo::Structure::BoundaryElementsCreate(int rElementGroupId, int rNodeGroupId, NodeBase* rControlNode)
{
    Timer timer(__FUNCTION__, GetShowTime(), GetLogger());

    // find groups
    boost::ptr_map<int, GroupBase>::iterator itGroupElements = mGroupMap.find(rElementGroupId);
    if (itGroupElements == mGroupMap.end())
        throw Exception(__PRETTY_FUNCTION__, "Group with the given identifier does not exist.");
    if (itGroupElements->second->GetType() != NuTo::eGroupId::Elements)
        throw Exception(__PRETTY_FUNCTION__, "Group is not an element group.");

    boost::ptr_map<int, GroupBase>::iterator itGroupBoundaryNodes = mGroupMap.find(rNodeGroupId);
    if (itGroupBoundaryNodes == mGroupMap.end())
        throw Exception(__PRETTY_FUNCTION__, "Group with the given identifier does not exist.");
    if (itGroupBoundaryNodes->second->GetType() != NuTo::eGroupId::Nodes)
        throw Exception(__PRETTY_FUNCTION__, "Group is not a node group.");

    Group<ElementBase>& elementGroup = *(itGroupElements->second->AsGroupElement());
    Group<NodeBase>& nodeGroup = *(itGroupBoundaryNodes->second->AsGroupNode());

    // since the search is done via the id's, the surface nodes are ptr, so make another set with the node ptrs
    std::set<const NodeBase*> nodePtrSet;
    for (auto itNode : nodeGroup)
    {
        nodePtrSet.insert(itNode.second);
    }

    std::vector<int> newBoundaryElementIds;

    // loop over all elements
    for (auto itElement : elementGroup)
    {
        ElementBase* elementPtr = itElement.second;
        const InterpolationType& interpolationType = elementPtr->GetInterpolationType();

        // std::cout << typeid(*elementPtr).name() << "\n";
        // std::cout << typeid(ContinuumElement<1>).name() << std::endl;
        if (typeid(*elementPtr) != typeid(ContinuumElement<1>) && typeid(*elementPtr) != typeid(ContinuumElement<2>) &&
            typeid(*elementPtr) != typeid(ContinuumElement<3>))
            throw Exception(__PRETTY_FUNCTION__, "Element is not a ContinuumElement.");

        // loop over all surfaces
        for (int iSurface = 0; iSurface < interpolationType.GetNumSurfaces(); ++iSurface)
        {
            bool elementSurfaceNodesMatchBoundaryNodes = true;
            Eigen::VectorXi surfaceNodeIndices = interpolationType.GetSurfaceNodeIndices(iSurface);

            int numSurfaceNodes = surfaceNodeIndices.rows();
            std::vector<const NodeBase*> surfaceNodes(numSurfaceNodes);

            for (int iSurfaceNode = 0; iSurfaceNode < numSurfaceNodes; ++iSurfaceNode)
            {
                surfaceNodes[iSurfaceNode] = elementPtr->GetNode(surfaceNodeIndices(iSurfaceNode, 0));
            }

            // check, if all surface nodes are in the node group
            for (auto& surfaceNode : surfaceNodes)
            {
                if (nodePtrSet.find(surfaceNode) == nodePtrSet.end())
                {
                    // this surface has at least one node that is not in the list, continue
                    elementSurfaceNodesMatchBoundaryNodes = false;
                    break;
                }
            }
            if (not elementSurfaceNodesMatchBoundaryNodes)
                continue;

            int surfaceId = iSurface;

            int elementId = GetUnusedId(mElementMap);

            ElementBase* boundaryElement = nullptr;
            ConstitutiveBase& constitutiveLaw = elementPtr->GetConstitutiveLaw(0);

            switch (elementPtr->GetLocalDimension())
            {
            case 1:
            {
                const auto& integrationType = *GetPtrIntegrationType(eIntegrationType::IntegrationType0DBoundary);
                auto& element = dynamic_cast<ContinuumElement<1>&>(*elementPtr);
                if (rControlNode == nullptr)
                    boundaryElement = new ContinuumBoundaryElement<1>(element, integrationType, surfaceId);
                else
                    boundaryElement = new ContinuumBoundaryElementConstrainedControlNode<1>(element, integrationType,
                                                                                            surfaceId, rControlNode);
                break;
            }
            case 2:
            {
                eIntegrationType integrationTypeEnum;
                // check for 2D types
                switch (interpolationType.GetStandardIntegrationType())
                {
                case eIntegrationType::IntegrationType2D3NGauss1Ip:
                case eIntegrationType::IntegrationType2D4NGauss1Ip:
                    integrationTypeEnum = eIntegrationType::IntegrationType1D2NGauss1Ip;
                    break;
                case eIntegrationType::IntegrationType2D3NGauss3Ip:
                case eIntegrationType::IntegrationType2D4NGauss4Ip:
//.........这里部分代码省略.........