Index: libraries/libmbmath/matvec6.h
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/matvec6.h,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -r1.16 -r1.17
--- libraries/libmbmath/matvec6.h	28 May 2007 21:03:20 -0000	1.16
+++ libraries/libmbmath/matvec6.h	28 Sep 2007 22:11:37 -0000	1.17
@@ -517,6 +517,20 @@
 		    m[0][0]*x.GetMat12()+m[0][1]*x.GetMat22(),
 		    m[1][0]*x.GetMat12()+m[1][1]*x.GetMat22());
    };  
+   
+   bool operator == (const Mat6x6& x) const {
+      return m[0][0] == x.GetMat11()
+		&& m[1][0] == x.GetMat21()
+		&& m[0][1] == x.GetMat12()
+		&& m[1][1] == x.GetMat22();
+   };
+   
+   bool operator != (const Mat6x6& x) const {
+      return m[0][0] != x.GetMat11()
+		|| m[1][0] != x.GetMat21()
+		|| m[0][1] != x.GetMat12()
+		|| m[1][1] != x.GetMat22();
+   };
    
    Mat6x6 Transpose(void) {
       return Mat6x6(m[0][0].Transpose(),
Index: mbdyn/base/constltp_impl.cc
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/constltp_impl.cc,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -r1.16 -r1.17
--- mbdyn/base/constltp_impl.cc	20 Jun 2007 19:44:27 -0000	1.16
+++ mbdyn/base/constltp_impl.cc	28 Sep 2007 20:52:44 -0000	1.17
@@ -114,6 +114,7 @@
 {
 	const char *s, *sOrig = HP.IsWord(CL1DWordSet);
 	if (sOrig == 0) {
+		/* default to linear elastic? */
 		s = "linear" "elastic";
 		sOrig = "";
 
@@ -137,8 +138,14 @@
 {
 	const char *s, *sOrig = HP.IsWord(CL3DWordSet);
 	if (sOrig == 0) {
+#if 0
 		s = "linear" "elastic";
 		sOrig = "";
+#else
+		silent_cerr("unknown constitutive law 3D type "
+			"at line " << HP.GetLineData() << std::endl);
+		throw DataManager::ErrGeneric();
+#endif
 
 	} else {
 		s = sOrig;
@@ -160,8 +167,14 @@
 {
 	const char *s, *sOrig = HP.IsWord(CL6DWordSet);
 	if (sOrig == 0) {
+#if 0
 		s = "linear" "elastic";
 		sOrig = "";
+#else
+		silent_cerr("unknown constitutive law 6D type "
+			"at line " << HP.GetLineData() << std::endl);
+		throw DataManager::ErrGeneric();
+#endif
 
 	} else {
 		s = sOrig;
Index: mbdyn/base/constltp_nlp.cc
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/constltp_nlp.cc,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- mbdyn/base/constltp_nlp.cc	15 Jun 2007 15:31:04 -0000	1.3
+++ mbdyn/base/constltp_nlp.cc	28 Sep 2007 22:13:38 -0000	1.4
@@ -48,6 +48,7 @@
 class NLPViscoElasticConstitutiveLaw
 : public ElasticConstitutiveLaw<T, Tder> {
 private:
+	ConstLawType::Type CLType;
 	Tder FDE0;
 	Tder FDEPrime0;
 	std::vector<const DifferentiableScalarFunction *> FDEsc;
@@ -56,11 +57,13 @@
 public:
 	NLPViscoElasticConstitutiveLaw(const TplDriveCaller<T>* pDC,
 		const T& PStress,
+		const ConstLawType::Type& cltype,
 		const Tder& fde0,
 		const std::vector<const DifferentiableScalarFunction *>& fdesc,
 		const Tder& fdeprime0,
 		const std::vector<const DifferentiableScalarFunction *>& fdeprimesc)
 	: ElasticConstitutiveLaw<T, Tder>(pDC, PStress),
+	CLType(cltype),
 	FDE0(fde0), FDEPrime0(fdeprime0),
 	FDEsc(fdesc), FDEPrimesc(fdeprimesc)
 	{
@@ -73,7 +76,7 @@
 	};
 
 	ConstLawType::Type GetConstLawType(void) const {
-		return ConstLawType::VISCOELASTIC;
+		return CLType;
 	};
 
 	virtual ConstitutiveLaw<T, Tder>* pCopy(void) const {
@@ -83,7 +86,7 @@
 		SAFENEWWITHCONSTRUCTOR(pCL, cl,
 			cl(ElasticConstitutiveLaw<T, Tder>::pGetDriveCaller()->pCopy(),
 				ElasticConstitutiveLaw<T, Tder>::PreStress,
-				FDE0, FDEsc, FDEPrime0, FDEPrimesc));
+				CLType, FDE0, FDEsc, FDEPrime0, FDEPrimesc));
 		return pCL;
 	};
 
@@ -129,10 +132,15 @@
 		T E = ElasticConstitutiveLaw<T, Tder>::Epsilon
 			- ElasticConstitutiveLaw<T, Tder>::Get();
 
-		ConstitutiveLaw<T, Tder>::F = ElasticConstitutiveLaw<T, Tder>::PreStress
-			+ FDE0*E + FDEPrime0*EpsPrime;
-		ConstitutiveLaw<T, Tder>::FDE = FDE0;
-		ConstitutiveLaw<T, Tder>::FDEPrime = FDEPrime0;
+		ConstitutiveLaw<T, Tder>::F = ElasticConstitutiveLaw<T, Tder>::PreStress;
+		if (CLType & ConstLawType::ELASTIC) {
+			ConstitutiveLaw<T, Tder>::F += FDE0*E;
+			ConstitutiveLaw<T, Tder>::FDE = FDE0;
+		}
+		if (CLType & ConstLawType::VISCOUS) {
+			ConstitutiveLaw<T, Tder>::F += FDEPrime0*EpsPrime;
+			ConstitutiveLaw<T, Tder>::FDEPrime = FDEPrime0;
+		}
 
 		for (unsigned i = 0; i < FDEsc.size(); i++) {
 
@@ -148,26 +156,28 @@
 			 */
 
 			doublereal dEps = E(i + 1);
-			doublereal dEpsPrime = EpsPrime(i + 1);
 
-			doublereal df1 = 0.;
-			doublereal df1DE = 0.;
-			doublereal df2 = 0.;
-			doublereal df2DE = 0.;
+			doublereal f = 0.;
+			doublereal fde = 0.;
 
-			if (FDEsc[i]) {
-				df1 = (*FDEsc[i])(dEps);
-				df1DE = FDEsc[i]->ComputeDiff(dEps);
+			if ((CLType & ConstLawType::ELASTIC) && FDEsc[i]) {
+				doublereal df1 = (*FDEsc[i])(dEps);
+				doublereal df1DE = FDEsc[i]->ComputeDiff(dEps);
+				f += df1*dEps;
+				fde += df1 + df1DE*dEps;
 			}
 
-			if (FDEPrimesc[i]) {
-				df2 = (*FDEPrimesc[i])(dEps);
-				df2DE = FDEPrimesc[i]->ComputeDiff(dEps);
+			if ((CLType & ConstLawType::VISCOUS) && FDEPrimesc[i]) {
+				doublereal df2 = (*FDEPrimesc[i])(dEps);
+				doublereal df2DE = FDEPrimesc[i]->ComputeDiff(dEps);
+				doublereal dEpsPrime = EpsPrime(i + 1);
+				f += df2*dEpsPrime;
+				fde += df2DE*dEpsPrime;
+				ConstitutiveLaw<T, Tder>::FDEPrime(i + 1, i + 1) += df2;
 			}
 
-			ConstitutiveLaw<T, Tder>::F(i + 1) += df1*dEps + df2*dEpsPrime;
-			ConstitutiveLaw<T, Tder>::FDE(i + 1, i + 1) += df1 + df1DE*dEps + df2DE*dEpsPrime;
-			ConstitutiveLaw<T, Tder>::FDEPrime(i + 1, i + 1) += df2;
+			ConstitutiveLaw<T, Tder>::F(i + 1) += f;
+			ConstitutiveLaw<T, Tder>::FDE(i + 1, i + 1) += fde;
 		}
 	};
 
@@ -180,6 +190,7 @@
 class NLPViscoElasticConstitutiveLaw<doublereal, doublereal>
 : public ElasticConstitutiveLaw<doublereal, doublereal> {
 private:
+	ConstLawType::Type CLType;
 	doublereal FDE0;
 	doublereal FDEPrime0;
 	std::vector<const DifferentiableScalarFunction *> FDEsc;
@@ -188,11 +199,13 @@
 public:
 	NLPViscoElasticConstitutiveLaw(const TplDriveCaller<doublereal>* pDC,
 		const doublereal& PStress,
+		const ConstLawType::Type& cltype,
 		const doublereal& fde0,
 		const std::vector<const DifferentiableScalarFunction *>& fdesc,
 		const doublereal& fdeprime0,
 		const std::vector<const DifferentiableScalarFunction *>& fdeprimesc)
 	: ElasticConstitutiveLaw<doublereal, doublereal>(pDC, PStress),
+	CLType(cltype),
 	FDE0(fde0), FDEPrime0(fdeprime0),
 	FDEsc(fdesc), FDEPrimesc(fdeprimesc)
 	{
@@ -205,7 +218,7 @@
 	};
 
 	ConstLawType::Type GetConstLawType(void) const {
-		return ConstLawType::VISCOELASTIC;
+		return CLType;
 	};
 
 	virtual ConstitutiveLaw<doublereal, doublereal>* pCopy(void) const {
@@ -215,7 +228,7 @@
 		SAFENEWWITHCONSTRUCTOR(pCL, cl,
 			cl(ElasticConstitutiveLaw<doublereal, doublereal>::pGetDriveCaller()->pCopy(),
 				ElasticConstitutiveLaw<doublereal, doublereal>::PreStress,
-				FDE0, FDEsc, FDEPrime0, FDEPrimesc));
+				CLType, FDE0, FDEsc, FDEPrime0, FDEPrimesc));
 		return pCL;
 	};
 
@@ -262,10 +275,15 @@
 			- ElasticConstitutiveLaw<doublereal, doublereal>::Get();
 		doublereal dEpsPrime = EpsPrime;
 
-		ConstitutiveLaw<doublereal, doublereal>::F = ElasticConstitutiveLaw<doublereal, doublereal>::PreStress
-			+ FDE0*dEps + FDEPrime0*dEpsPrime;
-		ConstitutiveLaw<doublereal, doublereal>::FDE = FDE0;
-		ConstitutiveLaw<doublereal, doublereal>::FDEPrime = FDEPrime0;
+		ConstitutiveLaw<doublereal, doublereal>::F = ElasticConstitutiveLaw<doublereal, doublereal>::PreStress;
+		if (CLType & ConstLawType::ELASTIC) {
+			ConstitutiveLaw<doublereal, doublereal>::F += FDE0*dEps;
+			ConstitutiveLaw<doublereal, doublereal>::FDE = FDE0;
+		}
+		if (CLType & ConstLawType::VISCOUS) {
+			ConstitutiveLaw<doublereal, doublereal>::F += FDEPrime0*dEpsPrime;
+			ConstitutiveLaw<doublereal, doublereal>::FDEPrime = FDEPrime0;
+		}
 
 		/*
 		 * each diagonal coefficient is:
@@ -278,24 +296,26 @@
 		 * f/epsPrime = f''
 		 */
 
-		doublereal df1 = 0.;
-		doublereal df1DE = 0.;
-		doublereal df2 = 0.;
-		doublereal df2DE = 0.;
+		doublereal f = 0.;
+		doublereal fde = 0.;
 
-		if (FDEsc[0]) {
-			df1 = (*FDEsc[0])(dEps);
-			df1DE = FDEsc[0]->ComputeDiff(dEps);
+		if ((CLType & ConstLawType::ELASTIC) && FDEsc[0]) {
+			doublereal df1 = (*FDEsc[0])(dEps);
+			doublereal df1DE = FDEsc[0]->ComputeDiff(dEps);
+			f += df1*dEps;
+			fde += df1 + df1DE*dEps;
 		}
 
-		if (FDEPrimesc[0]) {
-			df2 = (*FDEPrimesc[0])(dEps);
-			df2DE = FDEPrimesc[0]->ComputeDiff(dEps);
+		if ((CLType & ConstLawType::VISCOUS)  && FDEPrimesc[0]) {
+			doublereal df2 = (*FDEPrimesc[0])(dEps);
+			doublereal df2DE = FDEPrimesc[0]->ComputeDiff(dEps);
+			f += df2*dEpsPrime;
+			fde += df2DE*dEpsPrime;
+			ConstitutiveLaw<doublereal, doublereal>::FDEPrime += df2;
 		}
 
-		ConstitutiveLaw<doublereal, doublereal>::F += df1*dEps + df2*dEpsPrime;
-		ConstitutiveLaw<doublereal, doublereal>::FDE += df1 + df1DE*dEps + df2DE*dEpsPrime;
-		ConstitutiveLaw<doublereal, doublereal>::FDEPrime += df2;
+		ConstitutiveLaw<doublereal, doublereal>::F += f;
+		ConstitutiveLaw<doublereal, doublereal>::FDE += fde;
 	};
 
 	virtual void IncrementalUpdate(const doublereal& DeltaEps, const doublereal& EpsPrime = 0.) {
@@ -315,8 +335,6 @@
 	{
 		ConstitutiveLaw<T, Tder>* pCL = 0;
 
-		CLType = ConstLawType::VISCOELASTIC;
-
 		unsigned dim;
 		if (typeid(T) == typeid(doublereal)) {
 			dim = 1;
@@ -339,6 +357,7 @@
 		Tder FDE0(0.);
 		FDE0 = HP.Get(FDE0);
 
+		bool bElastic(FDE0 != Tder(0.));
 		std::vector<const DifferentiableScalarFunction *> FDEsc(dim);
 		for (unsigned i = 0; i < dim; i++) {
 			if (HP.IsKeyWord("null")) {
@@ -354,6 +373,7 @@
 						"must be differentiable" << std::endl);
 					throw ErrGeneric();
 				}
+				bElastic = true;
 			}
 		}
 
@@ -361,6 +381,7 @@
 		Tder FDEPrime0(0.);
 		FDEPrime0 = HP.Get(FDEPrime0);
 
+		bool bViscous(FDEPrime0 != Tder(0.));
 		std::vector<const DifferentiableScalarFunction *> FDEPrimesc(dim);
 		for (unsigned i = 0; i < dim; i++) {
 			if (HP.IsKeyWord("null")) {
@@ -376,6 +397,7 @@
 						"must be differentiable" << std::endl);
 					throw ErrGeneric();
 				}
+				bViscous = true;
 			}
 		}
 
@@ -385,9 +407,21 @@
 		T PreStrain(0.);
 		TplDriveCaller<T>* pTplDC = GetPreStrain(pDM, HP, PreStrain);
 
+		if (bElastic && bViscous) {
+			CLType = ConstLawType::VISCOELASTIC;
+		} else if (bElastic) {
+			CLType = ConstLawType::ELASTIC;
+		} else if (bViscous) {
+			CLType = ConstLawType::VISCOUS;
+		} else {
+			/* needs to be at least elastic... */
+			CLType = ConstLawType::ELASTIC;
+		}
+
 		typedef NLPViscoElasticConstitutiveLaw<T, Tder> L;
 		SAFENEWWITHCONSTRUCTOR(pCL, L,
 			L(pTplDC, PreStress,
+				CLType,
 				FDE0, FDEsc,
 				FDEPrime0, FDEPrimesc));
 
Index: mbdyn/base/constltp_nlsf.cc
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/constltp_nlsf.cc,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -r1.3 -r1.4
--- mbdyn/base/constltp_nlsf.cc	15 Jun 2007 15:31:04 -0000	1.3
+++ mbdyn/base/constltp_nlsf.cc	28 Sep 2007 22:13:38 -0000	1.4
@@ -48,6 +48,7 @@
 class NLSFViscoElasticConstitutiveLaw
 : public ElasticConstitutiveLaw<T, Tder> {
 private:
+	ConstLawType::Type CLType;
 	Tder FDE0;
 	Tder FDEPrime0;
 	std::vector<const DifferentiableScalarFunction *> FDEsc;
@@ -56,11 +57,13 @@
 public:
 	NLSFViscoElasticConstitutiveLaw(const TplDriveCaller<T>* pDC,
 		const T& PStress,
+		const ConstLawType::Type &cltype,
 		const Tder& fde0,
 		const std::vector<const DifferentiableScalarFunction *>& fdesc,
 		const Tder& fdeprime0,
 		const std::vector<const DifferentiableScalarFunction *>& fdeprimesc)
 	: ElasticConstitutiveLaw<T, Tder>(pDC, PStress),
+	CLType(cltype),
 	FDE0(fde0), FDEPrime0(fdeprime0),
 	FDEsc(fdesc), FDEPrimesc(fdeprimesc)
 	{
@@ -73,7 +76,7 @@
 	};
 
 	ConstLawType::Type GetConstLawType(void) const {
-		return ConstLawType::VISCOELASTIC;
+		return CLType;
 	};
 
 	virtual ConstitutiveLaw<T, Tder>* pCopy(void) const {
@@ -83,7 +86,7 @@
 		SAFENEWWITHCONSTRUCTOR(pCL, cl,
 			cl(ElasticConstitutiveLaw<T, Tder>::pGetDriveCaller()->pCopy(),
 				ElasticConstitutiveLaw<T, Tder>::PreStress,
-				FDE0, FDEsc, FDEPrime0, FDEPrimesc));
+				CLType, FDE0, FDEsc, FDEPrime0, FDEPrimesc));
 		return pCL;
 	};
 
@@ -129,10 +132,15 @@
 		T E = ElasticConstitutiveLaw<T, Tder>::Epsilon
 			- ElasticConstitutiveLaw<T, Tder>::Get();
 
-		ConstitutiveLaw<T, Tder>::F = ElasticConstitutiveLaw<T, Tder>::PreStress
-			+ FDE0*E + FDEPrime0*EpsPrime;
-		ConstitutiveLaw<T, Tder>::FDE = FDE0;
-		ConstitutiveLaw<T, Tder>::FDEPrime = FDEPrime0;
+		ConstitutiveLaw<T, Tder>::F = ElasticConstitutiveLaw<T, Tder>::PreStress;
+		if (CLType & ConstLawType::ELASTIC) {
+			ConstitutiveLaw<T, Tder>::F += FDE0*E;
+			ConstitutiveLaw<T, Tder>::FDE = FDE0;
+		}
+		if (CLType & ConstLawType::VISCOUS) {
+			ConstitutiveLaw<T, Tder>::F += FDEPrime0*EpsPrime;
+			ConstitutiveLaw<T, Tder>::FDEPrime = FDEPrime0;
+		}
 
 		for (unsigned i = 0; i < FDEsc.size(); i++) {
 
@@ -147,27 +155,23 @@
 			 * f/epsPrime = f0'' + f''/epsPrime
 			 */
 
-			doublereal dEps = E(i + 1);
-			doublereal dEpsPrime = EpsPrime(i + 1);
+			doublereal f = 0.;
 
-			doublereal df1 = 0.;
-			doublereal df1DE = 0.;
-			doublereal df2 = 0.;
-			doublereal df2DEPrime = 0.;
-
-			if (FDEsc[i]) {
-				df1 = (*FDEsc[i])(dEps);
-				df1DE = FDEsc[i]->ComputeDiff(dEps);
+			if ((CLType & ConstLawType::ELASTIC) && FDEsc[i]) {
+				doublereal dEps = E(i + 1);
+				f += (*FDEsc[i])(dEps);
+				doublereal df1DE = FDEsc[i]->ComputeDiff(dEps);
+				ConstitutiveLaw<T, Tder>::FDE(i + 1, i + 1) += df1DE;
 			}
 
-			if (FDEPrimesc[i]) {
-				df2 = (*FDEPrimesc[i])(dEpsPrime);
-				df2DEPrime = FDEPrimesc[i]->ComputeDiff(dEpsPrime);
+			if ((CLType & ConstLawType::ELASTIC) && FDEPrimesc[i]) {
+				doublereal dEpsPrime = EpsPrime(i + 1);
+				f += (*FDEPrimesc[i])(dEpsPrime);
+				doublereal df2DEPrime = FDEPrimesc[i]->ComputeDiff(dEpsPrime);
+				ConstitutiveLaw<T, Tder>::FDEPrime(i + 1, i + 1) += df2DEPrime;
 			}
 
-			ConstitutiveLaw<T, Tder>::F(i + 1) += df1 + df2;
-			ConstitutiveLaw<T, Tder>::FDE(i + 1, i + 1) += df1DE;
-			ConstitutiveLaw<T, Tder>::FDEPrime(i + 1, i + 1) += df2DEPrime;
+			ConstitutiveLaw<T, Tder>::F(i + 1) += f;
 		}
 	};
 
@@ -180,6 +184,7 @@
 class NLSFViscoElasticConstitutiveLaw<doublereal, doublereal>
 : public ElasticConstitutiveLaw<doublereal, doublereal> {
 private:
+	ConstLawType::Type CLType;
 	doublereal FDE0;
 	doublereal FDEPrime0;
 	std::vector<const DifferentiableScalarFunction *> FDEsc;
@@ -188,11 +193,13 @@
 public:
 	NLSFViscoElasticConstitutiveLaw(const TplDriveCaller<doublereal>* pDC,
 		const doublereal& PStress,
+		const ConstLawType::Type& cltype,
 		const doublereal& fde0,
 		const std::vector<const DifferentiableScalarFunction *>& fdesc,
 		const doublereal& fdeprime0,
 		const std::vector<const DifferentiableScalarFunction *>& fdeprimesc)
 	: ElasticConstitutiveLaw<doublereal, doublereal>(pDC, PStress),
+	CLType(cltype),
 	FDE0(fde0), FDEPrime0(fdeprime0),
 	FDEsc(fdesc), FDEPrimesc(fdeprimesc)
 	{
@@ -205,7 +212,7 @@
 	};
 
 	ConstLawType::Type GetConstLawType(void) const {
-		return ConstLawType::VISCOELASTIC;
+		return CLType;
 	};
 
 	virtual ConstitutiveLaw<doublereal, doublereal>* pCopy(void) const {
@@ -215,7 +222,7 @@
 		SAFENEWWITHCONSTRUCTOR(pCL, cl,
 			cl(ElasticConstitutiveLaw<doublereal, doublereal>::pGetDriveCaller()->pCopy(),
 				ElasticConstitutiveLaw<doublereal, doublereal>::PreStress,
-				FDE0, FDEsc, FDEPrime0, FDEPrimesc));
+				CLType, FDE0, FDEsc, FDEPrime0, FDEPrimesc));
 		return pCL;
 	};
 
@@ -262,10 +269,15 @@
 			- ElasticConstitutiveLaw<doublereal, doublereal>::Get();
 		doublereal dEpsPrime = EpsPrime;
 
-		ConstitutiveLaw<doublereal, doublereal>::F = ElasticConstitutiveLaw<doublereal, doublereal>::PreStress
-			+ FDE0*dEps + FDEPrime0*dEpsPrime;
-		ConstitutiveLaw<doublereal, doublereal>::FDE = FDE0;
-		ConstitutiveLaw<doublereal, doublereal>::FDEPrime = FDEPrime0;
+		ConstitutiveLaw<doublereal, doublereal>::F = ElasticConstitutiveLaw<doublereal, doublereal>::PreStress;
+		if (CLType & ConstLawType::ELASTIC) {
+			ConstitutiveLaw<doublereal, doublereal>::F += FDE0*dEps;
+			ConstitutiveLaw<doublereal, doublereal>::FDE = FDE0;
+		}
+		if (CLType & ConstLawType::VISCOUS) {
+			ConstitutiveLaw<doublereal, doublereal>::F += FDEPrime0*dEpsPrime;
+			ConstitutiveLaw<doublereal, doublereal>::FDEPrime = FDEPrime0;
+		}
 
 		/*
 		 * each diagonal coefficient is:
@@ -278,24 +290,21 @@
 		 * f/epsPrime = f''
 		 */
 
-		doublereal df1 = 0.;
-		doublereal df1DE = 0.;
-		doublereal df2 = 0.;
-		doublereal df2DEPrime = 0.;
+		doublereal f = 0.;
 
-		if (FDEsc[0]) {
-			df1 = (*FDEsc[0])(dEps);
-			df1DE = FDEsc[0]->ComputeDiff(dEps);
+		if ((CLType & ConstLawType::ELASTIC) && FDEsc[0]) {
+			f += (*FDEsc[0])(dEps);
+			doublereal df1DE = FDEsc[0]->ComputeDiff(dEps);
+			ConstitutiveLaw<doublereal, doublereal>::FDE += df1DE;
 		}
 
-		if (FDEPrimesc[0]) {
-			df2 = (*FDEPrimesc[0])(dEpsPrime);
-			df2DEPrime = FDEPrimesc[0]->ComputeDiff(dEpsPrime);
+		if ((CLType & ConstLawType::VISCOUS) && FDEPrimesc[0]) {
+			f += (*FDEPrimesc[0])(dEpsPrime);
+			doublereal df2DEPrime = FDEPrimesc[0]->ComputeDiff(dEpsPrime);
+			ConstitutiveLaw<doublereal, doublereal>::FDEPrime += df2DEPrime;
 		}
 
-		ConstitutiveLaw<doublereal, doublereal>::F += df1 + df2;
-		ConstitutiveLaw<doublereal, doublereal>::FDE += df1DE;
-		ConstitutiveLaw<doublereal, doublereal>::FDEPrime += df2DEPrime;
+		ConstitutiveLaw<doublereal, doublereal>::F += f;
 	};
 
 	virtual void IncrementalUpdate(const doublereal& DeltaEps, const doublereal& EpsPrime = 0.) {
@@ -315,8 +324,6 @@
 	{
 		ConstitutiveLaw<T, Tder>* pCL = 0;
 
-		CLType = ConstLawType::VISCOELASTIC;
-
 		unsigned dim;
 		if (typeid(T) == typeid(doublereal)) {
 			dim = 1;
@@ -339,6 +346,7 @@
 		Tder FDE0(0.);
 		FDE0 = HP.Get(FDE0);
 
+		bool bElastic(FDE0 != Tder(0.));
 		std::vector<const DifferentiableScalarFunction *> FDEsc(dim);
 		for (unsigned i = 0; i < dim; i++) {
 			if (HP.IsKeyWord("null")) {
@@ -354,6 +362,7 @@
 						"must be differentiable" << std::endl);
 					throw ErrGeneric();
 				}
+				bElastic = true;
 			}
 		}
 
@@ -361,6 +370,7 @@
 		Tder FDEPrime0(0.);
 		FDEPrime0 = HP.Get(FDEPrime0);
 
+		bool bViscous(FDEPrime0 != Tder(0.));
 		std::vector<const DifferentiableScalarFunction *> FDEPrimesc(dim);
 		for (unsigned i = 0; i < dim; i++) {
 			if (HP.IsKeyWord("null")) {
@@ -376,6 +386,7 @@
 						"must be differentiable" << std::endl);
 					throw ErrGeneric();
 				}
+				bViscous = true;
 			}
 		}
 
@@ -385,9 +396,21 @@
 		T PreStrain(0.);
 		TplDriveCaller<T>* pTplDC = GetPreStrain(pDM, HP, PreStrain);
 
+		if (bElastic && bViscous) {
+			CLType = ConstLawType::VISCOELASTIC;
+		} else if (bElastic) {
+			CLType = ConstLawType::ELASTIC;
+		} else if (bViscous) {
+			CLType = ConstLawType::VISCOUS;
+		} else {
+			/* needs to be at least elastic... */
+			CLType = ConstLawType::ELASTIC;
+		}
+
 		typedef NLSFViscoElasticConstitutiveLaw<T, Tder> L;
 		SAFENEWWITHCONSTRUCTOR(pCL, L,
 			L(pTplDC, PreStress,
+				CLType,
 				FDE0, FDEsc,
 				FDEPrime0, FDEPrimesc));
 
Index: mbdyn/base/mbpar.cc
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/mbpar.cc,v
retrieving revision 1.57
retrieving revision 1.58
diff -u -r1.57 -r1.58
--- mbdyn/base/mbpar.cc	13 Jun 2007 21:20:19 -0000	1.57
+++ mbdyn/base/mbpar.cc	28 Sep 2007 20:52:44 -0000	1.58
@@ -419,7 +419,7 @@
 		/* allow "reference" (copy cached constitutive law) */
 		ConstitutiveLaw3D *pCL = GetConstLaw3D(clt);
 		if (pCL == NULL) {
-			silent_cerr("unable to read constitutive law 1D " 
+			silent_cerr("unable to read constitutive law 3D " 
 					<< uLabel);
 			if (sName) {
 				silent_cerr(" (" << sName << ")");
@@ -1203,7 +1203,7 @@
 	unsigned int uLabel = GetInt();
 	C6DType::const_iterator i = C6D.find(uLabel);
 	if (i == C6D.end()) {
-		silent_cerr("constitutive law 3D " << uLabel
+		silent_cerr("constitutive law 6D " << uLabel
 				<< " is undefined at line "
 				<< GetLineData() << std::endl);
 		throw MBDynParser::ErrGeneric();