Now storing parts lists defined in input file as domain lists on the …

…mesh. Added new surface reconstruction algorithm for surfaces defined via part lists. Fixed some bugs in FEMeshTopo class.

FEAMR/FEErosionAdaptor.cpp

@@ -46,7 +46,7 @@ BEGIN_FECORE_CLASS(FEErosionAdaptor, FEMeshAdaptor)
 	ADD_PARAMETER(m_maxelem, "max_elems");
 	ADD_PARAMETER(m_nsort, "sort");
 	ADD_PARAMETER(m_bremoveIslands, "remove_islands");
-	ADD_PARAMETER(m_erodeSurfaces, "erode_surfaces")->setEnums("No\0Yes\0Grow\0");
+	ADD_PARAMETER(m_erodeSurfaces, "erode_surfaces")->setEnums("no\0yes\0grow\0reconstruct\0");
 
 	ADD_PROPERTY(m_criterion, "criterion");
 END_FECORE_CLASS();
@@ -136,6 +136,7 @@ bool FEErosionAdaptor::Apply(int iteration)
 	case SurfaceErodeOption::DONT_ERODE: break;// don't do anything
 	case SurfaceErodeOption::ERODE: ErodeSurfaces(); break;
 	case SurfaceErodeOption::GROW : GrowErodedSurfaces(topo); break;
+	case SurfaceErodeOption::RECONSTRUCT: ReconstructSurfaces(topo); break;
 	}
 
 	// remove any linear constraints of excluded nodes
@@ -178,7 +179,7 @@ void FEErosionAdaptor::RemoveIslands(FEMeshTopo& topo)
 				island.push_back(id);
 
 				// loop over all the neighbors
-				vector<int> nbrList = topo.ElementNeighborIndexList(n);
+				vector<int> nbrList = topo.ElementNeighborIndexList(id);
 				for (int i = 0; i < nbrList.size(); ++i)
 				{
 					FEElement* eli = topo.Element(nbrList[i]);
@@ -435,3 +436,40 @@ void FEErosionAdaptor::UpdateLinearConstraints()
 		if (del) LCM.RemoveLinearConstraint(j); else ++j;
 	}
 }
+
+void FEErosionAdaptor::ReconstructSurfaces(FEMeshTopo& topo)
+{
+	FEModel& fem = *GetFEModel();
+	FEMesh& mesh = fem.GetMesh();
+
+	// We assume that the surfaces are definded through part lists. 
+	// We don't actually store which surfaces are generated from part lists,
+	// but they should have the same name. 
+	for (int i = 0; i < mesh.Surfaces(); ++i)
+	{
+		FESurface& surf = mesh.Surface(i);
+
+		// see if a part list exists with this name
+		FEDomainList* domList = mesh.FindDomainList(surf.GetName());
+		if (domList)
+		{
+			// build the outer surface
+			FEFacetSet* facetSet = mesh.DomainBoundary(*domList);
+
+			FEFacetSet* oldFacetSet = surf.GetFacetSet();
+			if (oldFacetSet == nullptr)
+			{
+				surf.Create(*facetSet);
+			}
+			else
+			{
+				oldFacetSet->Clear();
+				oldFacetSet->Add(facetSet);
+				surf.Create(*oldFacetSet);
+			}
+
+			surf.CreateMaterialPointData();
+			surf.Init();
+		}
+	}
+}

FEAMR/FEErosionAdaptor.h

@@ -34,7 +34,8 @@ class FEErosionAdaptor : public FEMeshAdaptor
 	enum SurfaceErodeOption {
 		DONT_ERODE,		// don't erode surface elements (default)
 		ERODE,			// erode surface elements
-		GROW			// add newly exposed surface elements
+		GROW,			// add newly exposed surface elements
+		RECONSTRUCT		// reconstruct surfaces
 	};
 
 public:
@@ -47,6 +48,7 @@ class FEErosionAdaptor : public FEMeshAdaptor
 	void DeactivateOrphanedNodes();
 	void ErodeSurfaces();
 	void GrowErodedSurfaces(FEMeshTopo& topo);
+	void ReconstructSurfaces(FEMeshTopo& topo);
 	void UpdateLinearConstraints();
 
 private:

FEBioXML/FEBModel.cpp

@@ -679,6 +679,22 @@ bool FEBModel::BuildPart(FEModel& fem, Part& part, bool buildDomains, const Tran
 		fesurfPair->SetSecondarySurface(surf2);
 	}
 
+	// create domain lists
+	for (int i = 0; i < part.PartLists(); ++i)
+	{
+		FEBModel::PartList* partList = part.GetPartList(i);
+
+		FEDomainList* domList = new FEDomainList();
+		domList->SetName(partList->Name());
+		for (int j = 0; j < partList->Parts(); ++j)
+		{
+			FEDomain* dom = mesh.FindDomain(partList->PartName(j)); assert(dom);
+			if (dom) domList->AddDomain(dom);
+		}
+
+		mesh.AddDomainList(domList);
+	}
+
 	// create discrete element sets
 	int DSets = part.DiscreteSets();
 	for (int i = 0; i < DSets; ++i)

FEBioXML/FEBModel.h

@@ -207,6 +207,9 @@ class FEBModel
 		void SetPartList(const std::vector<std::string>& parts);
 		const std::vector<std::string>& GetPartList() const;
 
+		size_t Parts() const { return m_parts.size(); }
+		const std::string& PartName(size_t n) const { return m_parts[n]; }
+
 	private:
 		std::string					m_name;
 		std::vector<std::string>	m_parts;

FECore/FEDomainList.cpp

@@ -44,6 +44,11 @@ FEDomainList::FEDomainList(FEDomainList& domList)
 	m_dom = domList.m_dom;
 }
 
+FEDomainList::FEDomainList(std::vector<FEDomain*> domList)
+{
+	m_dom = domList;
+}
+
 //! Clear the domain list
 void FEDomainList::Clear()
 {

FECore/FEDomainList.h

@@ -23,11 +23,9 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.*/
-
-
-
 #pragma once
 #include <vector>
+#include <string>
 #include "fecore_api.h"
 
 class FEDomain;
@@ -41,6 +39,7 @@ class FECORE_API FEDomainList
 	//! constructor
 	FEDomainList();
 	FEDomainList(FEDomainList& domList);
+	FEDomainList(std::vector<FEDomain*> domList);
 
 	//! Clear the domain list
 	void Clear();
@@ -67,6 +66,14 @@ class FECORE_API FEDomainList
 	//! serialization
 	void Serialize(DumpStream& ar);
 
+	size_t size() const { return m_dom.size(); }
+	FEDomain* operator [] (size_t n) { return m_dom[n]; }
+
+public:
+	void SetName(const std::string& name) { m_name = name; }
+	const std::string& GetName() const { return m_name; }
+
 private:
+	std::string m_name;
 	std::vector<FEDomain*>	m_dom;		// the actual list of domains
 };

FECore/FEFacetSet.cpp

@@ -49,6 +49,11 @@ FEFacetSet::FEFacetSet(FEModel* fem) : FEItemList(fem, FEItemList::FE_ELEMENT_SE
 void FEFacetSet::SetSurface(FESurface* surf) { m_surface = surf; }
 FESurface* FEFacetSet::GetSurface() { return m_surface; }
 
+void FEFacetSet::Clear()
+{
+	m_Face.clear();
+}
+
 //-----------------------------------------------------------------------------
 int FEFacetSet::Faces() const { return (int)m_Face.size(); }
 

FECore/FEFacetSet.h

@@ -79,6 +79,8 @@ class FECORE_API FEFacetSet : public FEItemList
 	// create from a surface
 	void Create(const FESurface& surf);
 
+	void Clear();
+
 	// return the size of the facet ste
 	int Faces() const;
 

FECore/FEMesh.cpp

@@ -393,6 +393,7 @@ void FEMesh::Clear()
 	for (size_t i=0; i<m_DiscSet.size (); ++i) delete m_DiscSet [i];
 	for (size_t i=0; i<m_FaceSet.size (); ++i) delete m_FaceSet [i];
 	for (size_t i=0; i<m_SurfPair.size(); ++i) delete m_SurfPair[i];
+	for (size_t i=0; i<m_DomList.size (); ++i) delete m_DomList [i];
 
 	m_Domain.clear();
 	m_Surf.clear();
@@ -402,6 +403,7 @@ void FEMesh::Clear()
 	m_DiscSet.clear();
 	m_FaceSet.clear();
 	m_SurfPair.clear();
+	m_DomList.clear();
 
 	m_NEL.Clear();
 	m_EEL.Clear();
@@ -560,11 +562,17 @@ FEFacetSet* FEMesh::FindFacetSet(const std::string& name)
 	return 0;
 }
 
-//-----------------------------------------------------------------------------
 FESurfacePair* FEMesh::FindSurfacePair(const std::string& name)
 {
-	for (size_t i = 0; i<m_SurfPair.size(); ++i) if (m_SurfPair[i]->GetName() == name) return m_SurfPair[i];
-	return 0;
+	for (FESurfacePair* sp : m_SurfPair) if (sp->GetName() == name) return sp;
+	return nullptr;
+}
+
+FEDomainList* FEMesh::FindDomainList(const std::string& name)
+{
+	for (FEDomainList* dom : m_DomList) 
+		if (dom->GetName() == name) return dom;
+	return nullptr;
 }
 
 //-----------------------------------------------------------------------------
@@ -896,6 +904,12 @@ FESurface* FEMesh::ElementBoundarySurface(std::vector<FEDomain*> domains, bool b
 }
 
 FEFacetSet* FEMesh::DomainBoundary(std::vector<FEDomain*> domains, bool boutside, bool binside)
+{
+	FEDomainList tmp(domains);
+	return DomainBoundary(tmp, boutside, binside);
+}
+
+FEFacetSet* FEMesh::DomainBoundary(FEDomainList& domains, bool boutside, bool binside)
 {
 	if ((boutside == false) && (binside == false)) return nullptr;
 
@@ -913,14 +927,22 @@ FEFacetSet* FEMesh::DomainBoundary(std::vector<FEDomain*> domains, bool boutside
 		for (int j = 0; j < domains[i]->Elements(); j++)
 		{
 			FEElement& el = domains[i]->ElementRef(j);
-			int index = FindElementIndexFromID(el.GetID());
-			int nf = el.Faces();
-			for (int k = 0; k < nf; ++k)
+			if (el.isActive())
 			{
-				FEElement* pen = EEL.Neighbor(index, k);
-				if ((pen == nullptr) && boutside) ++NF;
-				else if (pen && (std::find(domains.begin(), domains.end(), pen->GetMeshPartition()) == domains.end()) && boutside) ++NF;
-				if ((pen != nullptr) && (el.GetID() < pen->GetID()) && binside && (std::find(domains.begin(), domains.end(), pen->GetMeshPartition()) != domains.end())) ++NF;
+				int index = FindElementIndexFromID(el.GetID());
+				int nf = el.Faces();
+				for (int k = 0; k < nf; ++k)
+				{
+					FEElement* pen = EEL.Neighbor(index, k);
+					if ((pen == nullptr) && boutside) ++NF;
+					else if (pen && !pen->isActive() && boutside) ++NF;
+					else if (pen)
+					{
+						FEDomain* domk = dynamic_cast<FEDomain*>(pen->GetMeshPartition()); assert(domk);
+						if (boutside && !domains.IsMember(domk)) ++NF;
+						else if (binside && domains.IsMember(domk) && (el.GetID() < pen->GetID())) ++NF;
+					}
+				}
 			}
 		}
 	}
@@ -938,33 +960,45 @@ FEFacetSet* FEMesh::DomainBoundary(std::vector<FEDomain*> domains, bool boutside
 		for (int j = 0; j < domains[i]->Elements(); j++)
 		{
 			FEElement& el = domains[i]->ElementRef(j);
-			int index = FindElementIndexFromID(el.GetID());
-			int nf = el.Faces();
-			for (int k = 0; k < nf; ++k)
+			if (el.isActive())
 			{
-				FEElement* pen = EEL.Neighbor(index, k);
-				if (((pen == nullptr) && boutside) ||
-					(pen && (std::find(domains.begin(), domains.end(), pen->GetMeshPartition()) == domains.end()) && boutside) ||
-					((pen != nullptr) && (el.GetID() < pen->GetID()) && binside && (std::find(domains.begin(), domains.end(), pen->GetMeshPartition()) != domains.end())))
+				int index = FindElementIndexFromID(el.GetID());
+				int nf = el.Faces();
+				for (int k = 0; k < nf; ++k)
 				{
-					FEFacetSet::FACET& f = ps->Face(NF++);
-					int fn = el.GetFace(k, faceNodes);
+					FEElement* pen = EEL.Neighbor(index, k);
+					bool addFace = false;
 
-					switch (fn)
+					if ((pen == nullptr) && boutside) addFace = true;
+					else if (pen && !pen->isActive() && boutside) addFace = true;
+					else if (pen)
 					{
-					case 4: f.ntype = FEFacetSet::FACET::QUAD4; break;
-					case 8: f.ntype = FEFacetSet::FACET::QUAD8; break;
-					case 9: f.ntype = FEFacetSet::FACET::QUAD9; break;
-					case 3: f.ntype = FEFacetSet::FACET::TRI3; break;
-					case 6: f.ntype = FEFacetSet::FACET::TRI6; break;
-					case 7: f.ntype = FEFacetSet::FACET::TRI7; break;
-					default:
-						assert(false);
+						FEDomain* domk = dynamic_cast<FEDomain*>(pen->GetMeshPartition()); assert(domk);
+						if (boutside && !domains.IsMember(domk)) addFace = true;
+						else if (binside && domains.IsMember(domk) && (el.GetID() < pen->GetID())) addFace = true;
 					}
 
-					for (int p = 0; p < fn; ++p)
+					if (addFace)
 					{
-						f.node[p] = faceNodes[p];
+						FEFacetSet::FACET& f = ps->Face(NF++);
+						int fn = el.GetFace(k, faceNodes);
+
+						switch (fn)
+						{
+						case 4: f.ntype = FEFacetSet::FACET::QUAD4; break;
+						case 8: f.ntype = FEFacetSet::FACET::QUAD8; break;
+						case 9: f.ntype = FEFacetSet::FACET::QUAD9; break;
+						case 3: f.ntype = FEFacetSet::FACET::TRI3; break;
+						case 6: f.ntype = FEFacetSet::FACET::TRI6; break;
+						case 7: f.ntype = FEFacetSet::FACET::TRI7; break;
+						default:
+							assert(false);
+						}
+
+						for (int p = 0; p < fn; ++p)
+						{
+							f.node[p] = faceNodes[p];
+						}
 					}
 				}
 			}

FECore/FEMesh.h

@@ -39,6 +39,7 @@ SOFTWARE.*/
 #include "FEBoundingBox.h"
 #include "FESolidElement.h"
 #include "FEShellElement.h"
+#include "FEDomainList.h"
 
 //-----------------------------------------------------------------------------
 class FEEdge;
@@ -245,6 +246,12 @@ class FECORE_API FEMesh
 	void AddSurfacePair(FESurfacePair* ps) { m_SurfPair.push_back(ps); }
 	FESurfacePair* FindSurfacePair(const std::string& name);
 
+public:
+	size_t DomainLists() const { return m_DomList.size(); }
+	FEDomainList* DomainList(size_t n) { return m_DomList[n]; }
+	void AddDomainList(FEDomainList* dl) { m_DomList.push_back(dl); }
+	FEDomainList* FindDomainList(const std::string& name);
+
 public:
 	//! stream mesh data
 	void Serialize(DumpStream& dmp);
@@ -267,6 +274,7 @@ class FECORE_API FEMesh
 	//! binside  : include all interior facets
 	FESurface* ElementBoundarySurface(std::vector<FEDomain*> domains, bool boutside = true, bool binside = false);
 	FEFacetSet* DomainBoundary(std::vector<FEDomain*> domains, bool boutside = true, bool binside = false);
+	FEFacetSet* DomainBoundary(FEDomainList& domainList, bool boutside = true, bool binside = false);
 
 	//! get the nodal coordinates in reference configuration
 	void GetInitialNodalCoordinates(const FEElement& el, vec3d* node);
@@ -300,6 +308,7 @@ class FECORE_API FEMesh
 	vector<FEDiscreteSet*>	m_DiscSet;	//!< discrete element sets
 	vector<FEFacetSet*>		m_FaceSet;	//!< facet sets
 	vector<FESurfacePair*>	m_SurfPair;	//!< facet set pairs
+	vector<FEDomainList*>	m_DomList;	//!< named domain lists
 
 	vector<FEDataMap*>		m_DataMap;	//!< all data maps
 

FECore/FEMeshTopo.cpp

@@ -290,7 +290,7 @@ std::vector<FEElement*> FEMeshTopo::ElementNeighborList(int n)
 	int nbrs = 0;
 	switch (el->Shape())
 	{
-	case ET_HEX8: nbrs = 8; break;
+	case ET_HEX8: nbrs = 6; break;
 	case ET_TET4: nbrs = 4; break;
 	case ET_TET5: nbrs = 4; break;
 	default:
@@ -313,7 +313,7 @@ std::vector<int> FEMeshTopo::ElementNeighborIndexList(int n)
 	int nbrs = 0;
 	switch (el->Shape())
 	{
-	case ET_HEX8: nbrs = 8; break;
+	case ET_HEX8: nbrs = 6; break;
 	case ET_TET4: nbrs = 4; break;
 	case ET_TET5: nbrs = 4; break;
 	default: