C++ Octree::elems方法代码示例

void Rendering::Mesh::init(FEMMesh* inputMesh)
{

	if (inputMesh)
	{
		LOGI("Creating mapping between simulation mesh \"%s\" and surface mesh" ,inputMesh->filename.c_str());
		static std::string input_filename;
		static sofa::helper::vector<Mat3x3d> bases;
		static sofa::helper::vector<Vec3d> centers;
		static Octree<Vec3d> octree;
		const TVecTetra& tetras = inputMesh->tetrahedra;
		const TVecCoord& in = inputMesh->positions;
		const TVecCoord& out = m_vertices;
		d_map_i.resize(out.size());
		d_map_f.resize(out.size());
		if (input_filename != inputMesh->filename || bases.size() != tetras.size()) // we have to recompute the octree and bases
		{
			input_filename = inputMesh->filename;
			sofa::helper::vector< BBox<Vec3d> > bbox;
			bases.resize(tetras.size());
			centers.resize(tetras.size());
			bbox.resize(tetras.size());
			LOGI("   Preparing tetrahedra");
			for (unsigned int t=0; t<tetras.size(); ++t)
			{
				Mat3x3d m, mt;
				m[0] = in[tetras[t][1]]-in[tetras[t][0]];
				m[1] = in[tetras[t][2]]-in[tetras[t][0]];
				m[2] = in[tetras[t][3]]-in[tetras[t][0]];
				mt.transpose(m);
				bases[t].invert(mt);
				centers[t] = (in[tetras[t][0]]+in[tetras[t][1]]+in[tetras[t][2]]+in[tetras[t][3]])*0.25;
				bbox[t].add(tetras[t].begin(), tetras[t].end(), in);
			}
			LOGI("  Building octree");
			octree.init(bbox,8,8);
		}
		LOGI( "  Processing vertices" );
		int outside = 0;
		sofa::helper::vector<Octree<Vec3d>*> cells;
		for (unsigned int i=0;i<out.size();i++)
		{
			Vec3d pos = out[i];
			Vec3d coefs;
			int index = -1;
			double distance = 1e10;
			Octree<Vec3d>* cell = octree.findNear(pos);
			if (cell)
			{
				const sofa::helper::vector<int>& elems = cell->elems();
				for (unsigned int e = 0; e < elems.size(); e++)
				{
					unsigned int t = elems[e];
					Vec3d v = bases[t] * (pos - in[tetras[t][0]]);
					double d = std::max(std::max(-v[0],-v[1]),std::max(-v[2],v[0]+v[1]+v[2]-1));
					if (d>0) d = (pos-centers[t]).norm2();
					if (d<distance) { coefs = v; distance = d; index = t; }
				}
			}
			if (distance > 0)
			{ // pos is outside of the fem mesh, find the nearest tetra

				// first let's find at least one tetra that is close, if not already found
				if (index >= 0) // we already have a close tetra, we need to look only for closer ones
				{
					cells.clear();
					octree.findAllAround(cells, pos, sqrt(distance)*1.5);
					for (unsigned int ci = 0; ci < cells.size(); ++ci)
					{
						if (cells[ci] == cell) continue; // already processed this cell
						const sofa::helper::vector<int>& elems = cells[ci]->elems();
						for (unsigned int e = 0; e < elems.size(); e++)
						{
							unsigned int t = elems[e];
							double d = (pos-centers[t]).norm2();
							if (d<distance)
							{
								coefs = bases[t] * (pos - in[tetras[t][0]]);
								distance = d; index = t;
							}
						}
					}
				}
				else
				{
					// failsafe case (should not happen...), to be sure we do a brute-force search
					for (unsigned int t = 0; t < tetras.size(); t++)
					{
						double d = (pos-centers[t]).norm2();
						if (d<distance)
						{
							coefs = bases[t] * (pos - in[tetras[t][0]]);
							distance = d; index = t;
						}
					}
				}
				if (index >= 0)
				{
					//if (verbose >= 1) std::cout << "Surface vertex " << i << " mapped outside of tetra " << index << " with coefs " << coefs << std::endl;
					++outside;
//.........这里部分代码省略.........

void SurfaceMesh::init(FEMMesh* inputMesh)
{
    // elements <-> particles table
#ifdef PARALLEL_GATHER
    {
        const int nbp = positions.size();
        const int nbe = triangles.size();
        const int nbBp = (nbp + BSIZE-1)/BSIZE;
        const int nbBe = (nbe + BSIZE-1)/BSIZE;
        // first find number of elements per particle
        std::vector<int> p_nbe;
        p_nbe.resize(nbp);
        for (int eindex = 0; eindex < nbe; ++eindex)
            for (int j = 0; j < 3; ++j)
                ++p_nbe[triangles[eindex][j]];
        // then compute max value
        nbElemPerVertex = 0;
        for (int i=0;i<nbp;++i)
            if (p_nbe[i] > nbElemPerVertex) nbElemPerVertex = p_nbe[i];
        // finally fill velems array
        velems.resize(nbBp*nbElemPerVertex*BSIZE);
        p_nbe.clear();
        p_nbe.resize(nbp);
        for (int eindex = 0; eindex < nbe; ++eindex)
            for (int j = 0; j < 3; ++j)
            {
                int p = triangles[eindex][j];
                int num = p_nbe[p]++;
                const int block  = p / BSIZE;
                const int thread = p % BSIZE;
                velems[ block * (nbElemPerVertex * BSIZE) +
                           num * BSIZE + thread ] = 1 + eindex;
            }
    }
#endif
	
    // FEMMesh -> SurfaceMesh mapping
    if (inputMesh)
    {
        std::cout << "Creating mapping between simulation mesh \"" << inputMesh->filename << "\" and surface mesh \"" << filename << "\"..." << std::endl;
        static std::string input_filename;
        static sofa::helper::vector<Mat3x3d> bases;
        static sofa::helper::vector<Vec3d> centers;
        static Octree<Vec3d> octree;
        const TVecTetra& tetras = inputMesh->tetrahedra;
        const TVecCoord& in = inputMesh->positions;
        const TVecCoord& out = positions;
        map_i.resize(out.size());
        map_f.resize(out.size());
        if (input_filename != inputMesh->filename || bases.size() != tetras.size()) // we have to recompute the octree and bases
        {
            input_filename = inputMesh->filename;
            sofa::helper::vector< BBox<Vec3d> > bbox;
            bases.resize(tetras.size());
            centers.resize(tetras.size());
            bbox.resize(tetras.size());
            std::cout << "  Preparing tetrahedra" << std::endl;
            for (unsigned int t=0; t<tetras.size(); ++t)
            {
                Mat3x3d m, mt;
                m[0] = in[tetras[t][1]]-in[tetras[t][0]];
                m[1] = in[tetras[t][2]]-in[tetras[t][0]];
                m[2] = in[tetras[t][3]]-in[tetras[t][0]];
                mt.transpose(m);
                bases[t].invert(mt);
                centers[t] = (in[tetras[t][0]]+in[tetras[t][1]]+in[tetras[t][2]]+in[tetras[t][3]])*0.25;
                bbox[t].add(tetras[t].begin(), tetras[t].end(), in);
            }
            std::cout << "  Building octree" << std::endl;
            octree.init(bbox,8,8);
        }
        std::cout << "  Processing vertices" << std::endl;
        int outside = 0;
        sofa::helper::vector<Octree<Vec3d>*> cells;
        for (unsigned int i=0;i<out.size();i++)
        {
            Vec3d pos = out[i];
            Vec3d coefs;
            int index = -1;
            double distance = 1e10;
            Octree<Vec3d>* cell = octree.findNear(pos);
            if (cell)
            {
                const sofa::helper::vector<int>& elems = cell->elems();
                for (unsigned int e = 0; e < elems.size(); e++)
                {
                    unsigned int t = elems[e];
                    Vec3d v = bases[t] * (pos - in[tetras[t][0]]);
                    double d = std::max(std::max(-v[0],-v[1]),std::max(-v[2],v[0]+v[1]+v[2]-1));
                    if (d>0) d = (pos-centers[t]).norm2();
                    if (d<distance) { coefs = v; distance = d; index = t; }
                }
            }
            if (distance > 0)
            { // pos is outside of the fem mesh, find the nearest tetra
                
                // first let's find at least one tetra that is close, if not already found
                if (index >= 0) // we already have a close tetra, we need to look only for closer ones
                {
                    cells.clear();
//.........这里部分代码省略.........