Fixed issues with incorrect use of material point classes.

FEBioFluid/FEBioFluidPlot.cpp

@@ -778,8 +778,8 @@ class FEFluidBodyForce
         for (int j = 0; j<NBL; ++j)
         {
             FEBodyForce* pbf = dynamic_cast<FEBodyForce*>(m_fem->ModelLoad(j));
-            FEMaterialPoint pt(mp);
-            if (pbf && pbf->IsActive()) bf += pbf->force(pt);
+			FEMaterialPoint& pt = const_cast<FEMaterialPoint&>(mp);
+			if (pbf && pbf->IsActive()) bf += pbf->force(pt);
         }
 		// FEBio actually applies the negative of the body force
         return -bf;

FEBioFluid/FEElasticFluid.cpp

@@ -229,7 +229,6 @@ double FEElasticFluid::Pressure(const double ef, const double T)
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double p = Pressure(tmp);
-    delete ft;
     return p;
 }
 
@@ -243,7 +242,6 @@ double FEElasticFluid::Tangent_Strain(const double ef, const double T)
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double dpJ = Tangent_Strain(tmp);
-    delete ft;
     return dpJ;
 }
 
@@ -257,7 +255,6 @@ double FEElasticFluid::Tangent_Temperature(const double ef, const double T)
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double dpT = Tangent_Temperature(tmp);
-    delete ft;
     return dpT;
 }
 
@@ -271,7 +268,6 @@ double FEElasticFluid::Tangent_Strain_Strain(const double ef, const double T)
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double dpJ2 = Tangent_Strain_Strain(tmp);
-    delete ft;
     return dpJ2;
 }
 
@@ -285,7 +281,6 @@ double FEElasticFluid::Tangent_Strain_Temperature(const double ef, const double
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double dpJT = Tangent_Strain_Temperature(tmp);
-    delete ft;
     return dpJT;
 }
 
@@ -299,6 +294,5 @@ double FEElasticFluid::Tangent_Temperature_Temperature(const double ef, const do
     ft->m_T = T;
     FEMaterialPoint tmp(ft);
     double dpT2 = Tangent_Temperature_Temperature(tmp);
-    delete ft;
     return dpT2;
 }

FEBioFluid/FEFluidSolutes.cpp

@@ -201,8 +201,8 @@ double FEFluidSolutes::ElectricPotential(const FEMaterialPoint& pt, const bool e
     const int nsol = (int)m_pSolute.size();
     double cF = 0.0;
 
-    FEMaterialPoint mp = pt;
-    vector<double> c(nsol);        // effective concentration
+	FEMaterialPoint& mp = const_cast<FEMaterialPoint&>(pt);
+	vector<double> c(nsol);        // effective concentration
     vector<double> khat(nsol);    // solubility
     vector<int> z(nsol);        // charge number
     for (i=0; i<nsol; ++i) {
@@ -248,7 +248,7 @@ double FEFluidSolutes::ElectricPotential(const FEMaterialPoint& pt, const bool e
 //! partition coefficient
 double FEFluidSolutes::PartitionCoefficient(const FEMaterialPoint& pt, const int sol)
 {
-    FEMaterialPoint mp = pt;
+    FEMaterialPoint& mp = const_cast<FEMaterialPoint&>(pt);
     // solubility
     double khat = m_pSolute[sol]->m_pSolub->Solubility(mp);
     // charge number
@@ -274,8 +274,8 @@ void FEFluidSolutes::PartitionCoefficientFunctions(const FEMaterialPoint& mp, ve
 
     const FEFluidMaterialPoint& fpt = *(mp.ExtractData<FEFluidMaterialPoint>());
     const FEFluidSolutesMaterialPoint& spt = *(mp.ExtractData<FEFluidSolutesMaterialPoint>());
-    FEMaterialPoint pt = mp;
-    
+	FEMaterialPoint& pt = const_cast<FEMaterialPoint&>(mp);
+
     const int nsol = (int)m_pSolute.size();
 
     vector<double> c(nsol);
@@ -435,7 +435,7 @@ double FEFluidSolutes::PressureActual(const FEMaterialPoint& pt)
         c[i] = ConcentrationActual(pt, i);
 
     // osmotic coefficient
-    FEMaterialPoint mp = pt;
+    FEMaterialPoint& mp = const_cast<FEMaterialPoint&>(pt);
     double osmc = m_pOsmC->OsmoticCoefficient(mp);
 
     // actual pressure
@@ -458,7 +458,7 @@ vec3d FEFluidSolutes::SoluteFlux(const FEMaterialPoint& pt, const int sol)
     vec3d gradc = spt.m_gradc[sol];
 
     // solute free diffusivity
-    FEMaterialPoint mp = pt;
+    FEMaterialPoint& mp = const_cast<FEMaterialPoint&>(pt);
     double D0 = m_pSolute[sol]->m_pDiff->Free_Diffusivity(mp);
     double kappa = PartitionCoefficient(pt, sol);
 
@@ -483,8 +483,8 @@ vec3d FEFluidSolutes::SoluteDiffusiveFlux(const FEMaterialPoint& pt, const int s
     vec3d gradc = spt.m_gradc[sol];
 
     // solute free diffusivity
-    FEMaterialPoint mp = pt;
-    double D0 = m_pSolute[sol]->m_pDiff->Free_Diffusivity(mp);
+	FEMaterialPoint& mp = const_cast<FEMaterialPoint&>(pt);
+	double D0 = m_pSolute[sol]->m_pDiff->Free_Diffusivity(mp);
     double kappa = PartitionCoefficient(pt, sol);
 
     // diffusive solute flux j

FEBioMech/FEBioMechPlot.cpp

@@ -1110,7 +1110,7 @@ class FESpecificStrainEnergy
 	double operator()(const FEMaterialPoint& mp)
 	{
         if (dynamic_cast<FERemodelingInterface*>(m_mat) == 0) return 0;
-        FEMaterialPoint lmp(mp);
+		FEMaterialPoint& lmp = const_cast<FEMaterialPoint&>(mp);
 		FERemodelingMaterialPoint* rpt = lmp.ExtractData<FERemodelingMaterialPoint>();
         if (rpt == nullptr) {
             FEMaterialPoint* pt = lmp.ExtractData<FEElasticMixtureMaterialPoint>()->GetPointData(m_comp);
@@ -4732,29 +4732,6 @@ bool FEPlotIdealGasPressure::Save(FESurface& surf, FEDataStream& a)
 	else return false;
 }
 
-class FEFluidBodyForce
-{
-public:
-	FEFluidBodyForce(FEModel* fem, FESolidDomain& dom) : m_fem(fem), m_dom(dom) {}
-
-	vec3d operator()(const FEMaterialPoint& mp)
-	{
-		int NBL = m_fem->ModelLoads();
-		vec3d bf(0, 0, 0);
-		for (int j = 0; j < NBL; ++j)
-		{
-			FEBodyForce* pbf = dynamic_cast<FEBodyForce*>(m_fem->ModelLoad(j));
-			FEMaterialPoint pt(mp);
-			if (pbf && pbf->IsActive()) bf += pbf->force(pt);
-		}
-		return bf;
-	}
-
-private:
-	FESolidDomain& m_dom;
-	FEModel* m_fem;
-};
-
 bool FEPlotBodyForce::Save(FEDomain& dom, FEDataStream& a)
 {
 	if (dom.Class() == FE_DOMAIN_SOLID)
@@ -4767,7 +4744,7 @@ bool FEPlotBodyForce::Save(FEDomain& dom, FEDataStream& a)
 			for (int j = 0; j < NBL; ++j)
 			{
 				FEBodyForce* pbf = dynamic_cast<FEBodyForce*>(fem->ModelLoad(j));
-				FEMaterialPoint pt(mp);
+				FEMaterialPoint& pt = const_cast<FEMaterialPoint&>(mp);
 				if (pbf && pbf->IsActive()) bf += pbf->force(pt);
 			}
 

FEBioMech/FEUT4Domain.cpp

@@ -361,8 +361,8 @@ void FEUT4Domain::Update(const FETimeInfo& tp)
 	//		 in other words, we loose the material axis orientation
 	//		 For now, I solve this by copying the Q parameter
 	//       from the first element that the node connects to
-	FEElasticMaterialPoint ep;
-	FEMaterialPoint mp(&ep);
+	FEElasticMaterialPoint* ep = new FEElasticMaterialPoint();
+	FEMaterialPoint mp(ep);
 	mp.Init();
 
 	// loop over all the nodes
@@ -374,8 +374,8 @@ void FEUT4Domain::Update(const FETimeInfo& tp)
 		mp.m_r0 = m_pMesh->Node(node.inode).m_r0;
 		mp.m_rt = m_pMesh->Node(node.inode).m_rt;
 
-		ep.m_F = node.Fi;
-		ep.m_J = ep.m_F.det();
+		ep->m_F = node.Fi;
+		ep->m_J = ep->m_F.det();
 
 		// calculate the stress
 		node.si = m_pMat->Stress(mp);
@@ -825,19 +825,19 @@ void FEUT4Domain::NodalMaterialStiffness(UT4NODE& node, matrix& ke, FESolidMater
 	int* peli = m_NEL.ElementIndexList(node.inode);
 
 	// create a material point
-	FEElasticMaterialPoint ep;
-	FEMaterialPoint mp(&ep);
+	FEElasticMaterialPoint* ep = new FEElasticMaterialPoint;
+	FEMaterialPoint mp(ep);
 	mp.Init();
 
 	// set the material point data
 	mp.m_r0 = m_pMesh->Node(node.inode).m_r0;
 	mp.m_rt = m_pMesh->Node(node.inode).m_rt;
 
-	ep.m_F = node.Fi;
-	ep.m_J = ep.m_F.det();
+	ep->m_F = node.Fi;
+	ep->m_J = ep->m_F.det();
 
 	// set the Cauchy-stress
-	ep.m_s = node.si;
+	ep->m_s = node.si;
 
 	// Calculate the spatial tangent
 	tens4ds C = pme->Tangent(mp);
@@ -847,15 +847,15 @@ void FEUT4Domain::NodalMaterialStiffness(UT4NODE& node, matrix& ke, FESolidMater
 	{
 		// subtract the isochoric component from C;
 		// C = C - a*Ciso = C - (a*(C - Cvol)) = (1-a)*C + a*Cvol
-		C = C*(1 - m_alpha) + Cvol(C, ep.m_s)*m_alpha;
+		C = C*(1 - m_alpha) + Cvol(C, ep->m_s)*m_alpha;
 	}
 	else
 	{
 		C = C*(1 - m_alpha);
 	}
 
 	// convert spatial matrix to material
-	C = spatial_to_material(C, ep.m_F);
+	C = spatial_to_material(C, ep->m_F);
 
 	// extract the 'D' matrix
 	double D[6][6] = {0};

FECore/FEMaterialPoint.cpp

@@ -92,6 +92,27 @@ FEMaterialPoint::~FEMaterialPoint()
 	delete m_data;
 }
 
+// NOTE: We need a copy constructor to create vectors of FEMaterialPoint,
+// but we should not actually need to copy anything
+FEMaterialPoint::FEMaterialPoint(const FEMaterialPoint&)
+{
+	m_elem = nullptr;
+	m_shape = nullptr;
+	m_data = nullptr;
+	m_J0 = m_Jt = 1;
+	m_index = -1;
+}
+
+FEMaterialPoint& FEMaterialPoint::operator = (const FEMaterialPoint&)
+{
+	m_elem = nullptr;
+	m_shape = nullptr;
+	m_data = nullptr;
+	m_J0 = m_Jt = 1;
+	m_index = -1;
+	return *this;
+}
+
 void FEMaterialPoint::Init()
 {
 	if (m_data) m_data->Init();

FECore/FEMaterialPoint.h

@@ -107,6 +107,11 @@ class FECORE_API FEMaterialPoint
 	FEMaterialPoint(FEMaterialPointData* data = nullptr);
 	virtual ~FEMaterialPoint();
 
+	// TODO: These functions copy  nothing! They are only included because we need them to create vectors!
+	//       I would like to delete these functions, but this means they cannot be used in vectors anymore.
+	FEMaterialPoint(const FEMaterialPoint&);
+	FEMaterialPoint& operator = (const FEMaterialPoint&);
+
 	//! The init function is used to intialize data
 	virtual void Init();