/*
 $Id: problemtr.cc,v 1.5 1996/11/20 10:00:10 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "problemtr.h"
#include "nodeco.h"

#include "physics.h"

#include "triang.h"
#include "elements.h"
#include "dirichlettr.h"
#include "materialstr.h"

#include "adapt.h"
#include "linsystem.h"
#include "sysmatml.h"
#include "sysmatbl.h"
#include "int.h"

#include "cmdpars.h"
extern CmdPars Cmd;

//-------------------------------------------------------------------------


TransientProblem:: TransientProblem() : uPrev(1), uPrevOnNewMesh(1)
{
    transientMode = True;

    lumpedMass = True; 		Cmd.get("lumpedMass", &lumpedMass);

    plotTimeStep = False;	Cmd.get("plotTimeStep", &plotTimeStep);
    postTimeStep = False;	Cmd.get("postTimeStep", &postTimeStep);
    writeTimeStep = False;	Cmd.get("writeTimeStep",&writeTimeStep);

    maxTimeSteps=1000;		Cmd.get("maxTimeSteps", &maxTimeSteps);

    startTime = 0.0;		Cmd.get("startTime",    &startTime);
    if (!Cmd.get("endTime",&endTime))  MissingParameter("endTime");


    theta = 1.0;		Cmd.get("theta",&theta);
    if (theta < 0.5 || theta > 1)
    {
	cout << "\n*** Wrong value for theta given (" << theta 
	<< ") -- must be between 0.5 and 1.0\n";
	cout.flush(); abort();
    }
    if (equal(theta,1.0)) implicitEuler = True;
    else 		  implicitEuler = False;

    midPointRule = False; 	Cmd.get("midPointRule", &midPointRule);
    if (midPointRule) { theta = 0.5; implicitEuler = False; }


    if ((constTimeStep = Cmd.isTrue("constTimeStep")))
    { 
	if (!Cmd.get("tau", &tau)) missingTau(); 
	if (tau >= endTime) endTime = tau;
    }

    if (!constTimeStep)
    {
	tau = 0.1*(endTime-startTime);  Cmd.get("Tau",&tau);  // init. guess

	staticFirstStep = Cmd.isTrue("staticFirstStep");

	if ((fixedFirstStep = Cmd.isTrue("fixedFirstStep")))
	{ 
	    if (!Cmd.get("tau", &tau)) missingTau(); 
	    if (tau >= endTime) tau = 0.99*endTime;
	}

	maxTau = 0.1*(endTime-startTime);   Cmd.get("maxTau", &maxTau);

	rhoTime  = 0.5;			Cmd.get("rhoTime",  &rhoTime); // 0.18
	rhoSpace = 1.0-rhoTime-0.1; 	Cmd.get("rhoSpace", &rhoSpace);
	fSpace   = 1.0;			Cmd.get("fSpace",   &fSpace);
	rhoSpace *= fSpace;

	maxReductions = 25;		Cmd.get("maxReductions", &maxReductions);

	dynamicScaling = Cmd.isSet("scaling","dynamic"); 
    }
}
//-------------------------------------------------------------------------

void TransientProblem:: missingTau()
{
    cout << "\n*** time step tau not specified\n";  cout.flush(); abort();
}
//-------------------------------------------------------------------------

void TransientProblem:: notImplemented(const char* name) const
{
      cout << "\n*** class TransientProblem: function " << name 
           << " not implemented\n";
      cout.flush(); abort();
}
//-------------------------------------------------------------------------

void TransientProblem:: assembleGlobalDefect(MLMatrix& /*AD*/, 
					     Vector<Num>& /*MD*/, Vector<Num>& /*bD*/)
{
    notImplemented("assembleGlobalDefect");
}
//-------------------------------------------------------------------------


Bool TransientProblem:: solve(Real relGlobalPrecision)
{ 
    if (constTimeStep) return fixedTimeSteps   (relGlobalPrecision);
    else	       return adaptiveTimeSteps(relGlobalPrecision);
}
//-------------------------------------------------------------------------


Bool TransientProblem:: fixedTimeSteps(Real relGlobalPrecision)
{ 
    Bool lastStep = False, absPrecision = False;	
    int  i, timeStep;

    time = startTime + tau;


    for (timeStep=1; timeStep<=maxTimeSteps; ++timeStep)
    {
	staticSolution(relGlobalPrecision, absPrecision, time);

	TransientSolutionInfo(timeStep, time, tau);

	if (lastStep) break;
	else if (time+tau >= endTime) { tau = endTime-time; lastStep = True; }
	if (equal(endTime+tau, endTime)) break;

	time += tau;

	uPrev.resize(u.high());
	FORALL(u,i) uPrev[i] = u[i];
	shiftMesh();
    }

    if (timeStep > maxTimeSteps)
    {
	cout << "\n** Time step limit reached: " << timeStep << "\n";
	return False;
    }
    else return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool TransientProblem:: adaptiveTimeSteps(Real globalPrecision)
{ 
    int  i, timeStep, reductions;
    Real statPrecision, newTau;
    Real UNorm, timeError, newUNorm;
    Bool lastStep = False;
    const Bool absPrecision = True;

    const Real TOL = sqrt(globalPrecision);

    time = startTime + tau; 	

    if 	    (staticFirstStep) staticInitialSolution(&UNorm, TOL);
    else if (fixedFirstStep)  initialUNorm(&UNorm);
    else		      adaptiveInitialSolution(&UNorm, TOL);

    newUNorm = UNorm;

    for (timeStep=1; timeStep<=maxTimeSteps; ++timeStep)
    {
	for (reductions=0; True; ++reductions)
	{
	    if (dynamicScaling)     UNorm = newUNorm;
	    else if (timeStep == 1) UNorm = Max(UNorm,newUNorm);

	    statPrecision = sqr(rhoSpace*TOL*UNorm);
	    staticSolution(statPrecision, absPrecision, time);

	    compTimeError(&timeError, &newUNorm, TOL); 

	    if (timeConvergence(TOL, timeError, UNorm, &newTau, timeStep,
				reductions))  break;

	    if (fixedFirstStep && timeStep==1) break;

	    time -= tau;
	    tau  =  newTau;
	    time += tau;
	}

	tau = newTau;
	TransientSolutionInfo(timeStep, time, tau);

	if (lastStep) break; 
	else if (time+tau >= endTime) { tau = endTime-time; lastStep = True; }
	if (equal(endTime+tau, endTime)) break;

	time += tau;

	uPrev.resize(u.high());
	FORALL(u,i) uPrev[i] = u[i];		 // save solution for next step
	shiftMesh();
    }

    if (timeStep > maxTimeSteps)
    {
	cout << "\n** Time step limit reached: " << timeStep << "\n";
	return False;
    }
    else return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// -- 	one coarse initial solution step to compute UNorm for scaling:


void TransientProblem:: initialUNorm(Real* UNorm)
{
    staticSolution(0.05, False, time);		// 5 % are sufficient
    
    *UNorm = compUNorm(u);

    cout << Form("\n Initial UNorm = %1.3g\n", *UNorm);
    TransientSolutionInfo(0, time, tau);
}
//-------------------------------------------------------------------------


// -- 	try to find a good guess for the first time step by successively 
//							reducing the tolerance


void TransientProblem:: adaptiveInitialSolution(Real* UNorm, Real TOL)
{
    Real tol, initTol, statPrecision, timeError, newTau;
    Bool absPrecision = True;
    int reductions;

    initialUNorm(UNorm);

    Real newUNorm = *UNorm;

    if (TOL >= 0.1) initTol = 2.1*TOL;
    else	    initTol = 0.1;

    cout << Form("\n Search for initial time step")
         << Form(" - start with TOL = %1.3g, UNorm = %1.3g\n",initTol,*UNorm);


    for (tol=initTol; tol >= 2.001*TOL; tol=0.5*tol)
    {
	for (reductions=0; True; ++reductions)
	{
	    *UNorm = Max(*UNorm,newUNorm); 
		
	    statPrecision = sqr(rhoSpace*tol*(*UNorm));
	    staticSolution(statPrecision, absPrecision, time);
		
	    compTimeError(&timeError, &newUNorm, tol); 
	    if (timeConvergence(tol, timeError, *UNorm, &newTau, 0, reductions))
									  break;
	    tau  = newTau;
	    time = startTime + tau;
	}
	    
	tau = newTau;
	    
	TransientSolutionInfo(0, time, tau);
	time = startTime + tau;
    }
}
//-------------------------------------------------------------------------


// -- 		static solution for tau = 0 to supply initial values


void TransientProblem:: staticInitialSolution(Real* UNorm, Real TOL)
{
    int  i;
    Bool absPrecision = False;
    Real precision;

    time = startTime; 	

    transientMode = False;

    // --  one first coarse solution to get UNorm for absolute precision

    staticSolution(0.05, absPrecision, time);
    *UNorm = compUNorm(u);

    cout << Form("\n Static Initial Solution: UNorm = %1.3g\n", *UNorm);

    absPrecision = True;
    precision = sqr(rhoSpace*TOL*(*UNorm));

    staticSolution(precision, absPrecision, time);
	
    TransientSolutionInfo(0, time, 0.0);

    time = startTime + tau;

    uPrev.resize(u.high());
    FORALL(u,i) uPrev[i] = u[i];		 // save solution for next step
    shiftMesh();
 
    transientMode = True;
}
//-------------------------------------------------------------------------


Real TransientProblem:: compUNorm(Vector<Num>& v, Vector<Num>* MDiag)
{
    int i;
    
    Real minU = ::real(v[1]);

    FORALL(v,i) { if (::real(v[i]) < minU)  minU = ::real(v[i]); }
    // FORALL(v,i) { if (::real(v[1]) < minU)  minU = ::real(v[1]); }

    FORALL(v,i) v[i] -= minU;
		  
    Real UNorm = 0.0;

    if (MDiag)  FORALL(v,i) UNorm += Abs((*MDiag)[i] * v[i]*v[i]);
    else  UNorm = 2.0*L2MassNorm(v);
		
    FORALL(v,i) v[i] += minU;
		  
    return sqrt(UNorm);
}
//-------------------------------------------------------------------------
Real TransientProblem:: L2MassNorm(const Vector<Num>& x) const 
{
    int  dim;
    Real L2 = 0.0;
    PATCH* patch;

    dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<int> globalNodes(dim); 
    Jacobian Jac;
    
    Mesh()->resetElemIter();

    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem);

	element->assembleP(*patch, AElem, Jac, 0);
	interface->getGlobalNodes(patch, globalNodes);  

	L2 += energyProd(AElem, x, globalNodes);
    }

    // L2 = sqrt(2.0*L2);		// annihilate factor 1/2 in energy
    return L2;
}
//-------------------------------------------------------------------------


Bool TransientProblem:: timeConvergence(Real TOL, Real timeError, Real UNorm, 
					Real* newTau, int timeStep,
					int reductions)
{
    Bool status;
    Real totalError, spaceError;

    if (timeError==0) *newTau = maxTau;
    else *newTau = tau * sqrt(rhoTime*TOL*UNorm/timeError); 

    if (*newTau > maxTau) *newTau = maxTau;

    spaceError = sqrt(2.0*Error.Top());
    totalError = timeError + spaceError;

  
    if ( (totalError < TOL*UNorm) || (totalError==0) ) status = True;  
    else
    {
	status = False;
	*newTau *= 0.5;

	if (reductions > maxReductions)
	{
	    cout << "\n*** Time Step reductions > maxReductions (" 
	         << maxReductions << ")\n"; cout.flush(); abort();
	}
    }

    timeErrorInfo(totalError, timeError, spaceError, UNorm, TOL, tau, *newTau, 
		  timeStep, reductions);

    return status;
}
//-------------------------------------------------------------------------


void TransientProblem:: timeErrorInfo(Real /*totalError*/, Real timeError, 
				      Real spaceError, Real UNorm, Real TOL, 
				      Real tau0, Real newTau, int timeStep,
				      int reductions)
{
    const char *format0 = "\n% 5s %10s %10s %7s %7s %8s %8s %3s %6s";
    const char *format  = 
            "\n %5d %8.3g %% %8.3g %% %7.3g %7.3g %8.3g %8.3g %3d %6d\n";

    Real absTOL = UNorm*TOL;

    cout << Form(format0, " TStep", "TimeError", "SpaceError", "TOL",
		   "Norm", "tau", "new/tau", "red", "nodes");

    cout << Form(format, timeStep, quot(100*timeError,absTOL),
		   quot(100*spaceError,absTOL), TOL, UNorm, tau0, 
		   quot(newTau,tau0), reductions, interface->Dim());

    cout << ' '; for (int i=1; i<=76; ++i) cout << '~';  cout << "\n";
    // Pause();
}
//-------------------------------------------------------------------------


void TransientProblem:: TransientSolutionInfo(int timeStep, Real time0, 
					      Real tau0)
{
    const int step = Error.high();

    Real error = sqrt(quot(Error.Top(),Energy.Top()));
    
    const char *format0 = "\n %8s %10s %8s %5s %5s %9s %10s %10s";
    const char *format  = "\n %8d %10.3g %8.3g %5d %5d %9d %10.6g %10.3g\n";

    cout << Form(format0,"TimeStep", "time(sec)", " new tau",
		   "Steps", "level", "nodes", "ENorm", "rel.Error");

    cout << Form(format, timeStep, time0, tau0, 
		   step, interface->MaxDepth(), interface->Dim(), 
		   sqrt(2.0*Energy.Top()), error);

    cout << ' '; for (int i=1; i<=76; ++i) cout << '=';  cout << "\n";

    if (plotTimeStep)  interface->plot (u, step, fileName, timeStep);
    if (writeTimeStep) interface->write(u, step, fileName);

    if (postTimeStep > 0)
      if (timeStep%postTimeStep == 0)
        interface->autoPost(u, step, fileName, timeStep);

    while ( Pause() == picture )
	interface->post (u, step, fileName);    

//    Pause();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// -- transport Solution of previous time step to the nodes of the new mesh


void TransientProblem:: solutionToNewMesh(Vector<Num>& uPrevOnNewMesh0, 
					  const Interface* interfaceDLY) const
{
    Timer timer;
    int   i, k, n, node, comp;
    Vector<Real> x(spaceDim);

    const Interface*  interf = interface;
    if (interfaceDLY) interf = interfaceDLY;

    const int lowNode  = uPrevOnNewMesh0.low();
    const int highNode = uPrevOnNewMesh0.high();

    NodeCoordinates  nc(spaceDim, lowNode, highNode);

    interf->getNodeCoordinates(nc);


    if (PrevMesh() == 0)			// the first time step
    {
	node = lowNode;
	for (n=lowNode; n<=highNode; n+=nComp)
	{
	    nc.getCoordinate(n,x);
	    
	    for (comp=1; comp<=nComp; ++comp)
	     uPrevOnNewMesh0[node++] = dirichletBCs->initialValue(x, startTime, comp);
	}
    }
    else
    {
	Vector<Real> xUnit(spaceDim);
	Vector<int>  prevNodes(element->NoOfNodes());
	Vector<Num>  uLocal(element->NoOfNodes());
	Vector<int>  nodes(interf->element->NoOfNodes());

	const PATCH* patch, *prevPatch;

	Vector<Bool> done(lowNode, highNode);
	FORALL(done,i) done[i] = False;


	MESH* mesh = Mesh();
	mesh->resetElemIter();

	while ((patch = mesh->elemIterAll()))
	{
	    interf->getGlobalNodes(patch, nodes);

	    for (n=1; n<=nodes.high(); n+=nComp)
	    {
		if ((node=nodes[n]) >= lowNode)
		{
		    if (!done[node])
		    {
			nc.getCoordinate(node, x);
			prevPatch = PrevMesh()->findPatch(x, xUnit, patch);
			
			interface->getGlobalNodes(prevPatch, prevNodes);

			FORALL(uLocal,k) uLocal[k] = uPrev[prevNodes[k]];
			
			if (nComp==1) uPrevOnNewMesh0[node] = 
						element->valueAt(xUnit, uLocal);
			else  element->valueAt(xUnit, uLocal, uPrevOnNewMesh0, 
					       node, nComp);
			done[node] = True;
		    }
		}
	    }
	}
    }

    if (Cmd.isTrue("TimeTransport")) 
      { cout << "\n\tSolution Transport: "; timer.cpu(); }

	/*
	FORALL(uPrevOnNewMesh0,node)
	{
	    nc.getCoordinate(node,x);
	    const PATCH* patch = PrevMesh()->findPatch(x, xUnit, 0);
	    
	    interface->getGlobalNodes(patch, prevNodes);
	    FORALL(prevNodes,i)  uLocal[i] = uPrev[prevNodes[i]];
	    
	    uPrevOnNewMesh0[node] = element->valueAt(xUnit, uLocal);
	}
	*/
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


TransientHeatConduction:: TransientHeatConduction() { }

//-------------------------------------------------------------------------
    

DirichletBCs* TransientHeatConduction:: newDirichletBCs()
{
    if (Cmd.isSet("DirichletBCs","ConstDirichlet"))
         return new ConstDirichletBCs(fileName);
    else if (Cmd.isSet("DirichletBCs","ConstTransDirichlet"))
         return new ConstTransDirichlet(fileName);
    else if (Cmd.isSet("DirichletBCs","StepDirichlet"))
         return new StepDirichlet();
    else if (Cmd.isSet("DirichletBCs","UserTransient"))
         return new UserTransDirichlet();
    else if (Cmd.isSet("DirichletBCs","JumpDirichlet"))
         return new JumpDirichlet();
    else if (Cmd.isSet("DirichletBCs","TransDirichlet"))
         return new TransDirichlet();
    else MissingParameter("DirichletBCs");
    return 0;
}
//-------------------------------------------------------------------------


void TransientHeatConduction:: newMaterial()
{
    delete material;

    if (Cmd.isSet("Material","UserTransient"))
    	 material = new UserTransMaterial(fileName, spaceDim);

    else if (Cmd.isSet("Material","StepSource"))
    	 material = new StepSource(fileName, spaceDim);

    else if (Cmd.isSet("Material","DefaultMaterial"))
    	 material = new DefaultMaterial(fileName, spaceDim);

    else if (Cmd.isSet("Material","TransPeakSource"))
    	 material = new TransPeakSource(fileName, spaceDim);

    else MissingParameter("Material");
}
//-------------------------------------------------------------------------


void TransientHeatConduction:: assembleGlobal(Real time0)
{
    if (transientMode)			// initial static solution possible
    {
	uPrevOnNewMesh.resize(1,interface->Dim());
	solutionToNewMesh(uPrevOnNewMesh);
    }
    Problem::assembleGlobal(time0);
}
//-------------------------------------------------------------------------


void TransientHeatConduction:: assemble(const Element& elem, const PATCH& t, 
					const Jacobian& Jac, Matrix<Num>& A, 
					Vector<Num>& b, 
					const Matrix<Bool>* pattern,
					Bool errorEstimatorCall)
{
    if (!transientMode) 		// initial static Solution
    {
	elem.assembleEllip (t,A,Jac,pattern);

	if (lumpedMass) elem.assembleLumpedMass(t,A,Jac,pattern);
	else  		elem.assembleMass(t,A,Jac,pattern);

	elem.assembleSource(t,b,Jac,pattern,time);
	
	if (t.onBoundary())
	{
	    elem.assembleCauchyBCs(t,A,b,pattern,time);
	    elem.assembleNeumannBCs(t,b,pattern,time);
	}

	//if (Cmd.isTrue("innerBoundary")) elem.assembleInnerBCs(t,b,pattern,time);
	if (Mesh()->innerBoundary)       elem.assembleInnerBCs(t,b,pattern,time);

	return;
    }

    int i, k;
    Real factor;

    const int dim = elem.NoOfNodes();
    Matrix<Real> E(dim,dim);		 // 'elliptic' sub-matrix
    Matrix<Real> P(dim,dim);		 // 'parabolic' sub-matrix (mass-matrix)
    Vector<Real> s(dim), s0(dim);


    for (i=1; i<=dim; ++i)
    for (k=1; k<=i;   ++k)  E(i,k) = P(i,k) = 0.0;


    if (!elem.assembleEllip(t,E,Jac,pattern)) 
      missingMaterialTerm("elliptic", t.Class());

    if (lumpedMass) elem.assembleLumpedMass(t,E,Jac,pattern);
    else  	    elem.assembleMass(t,E,Jac,pattern);

    if (!errorEstimatorCall) 
    {
	Bool pTerm = False;

	if (lumpedMass) pTerm = elem.assembleLumpedP(t,P,Jac,pattern);
	else  		pTerm = elem.assembleP(t,P,Jac,pattern);

	if (!pTerm) missingMaterialTerm("mass (heat capacitance)",t.Class());
    }
    else elem.assembleP(t,P,Jac,pattern);      // no lumping for error estimator!


    // -- the contribution of the previous step to the right-hand side:


    if (!errorEstimatorCall) 			// ! A is used as working space
    {
	Vector<int> nodes(dim);
	interface->getGlobalNodes(&t,nodes);	

	if (implicitEuler) 
	{ 
	    if (lumpedMass) FORALL(b,i) b[i] += P(i,i)*uPrevOnNewMesh[nodes[i]]; 
	    else
	    {
		for (i=1; i<=dim; ++i)
		for (k=1; k< i;   ++k) P(k,i) = P(i,k);

		for (i=1; i<=dim; ++i)
		for (k=1; k<=dim; ++k) b[i] += P(i,k)*uPrevOnNewMesh[nodes[k]]; 
	    }
	}
	else        
	{ 
	    factor = tau*(1.0-theta);

	    for (i=1; i<=dim; ++i)
	    for (k=1; k<=i;   ++k)  
	    { 
		A(i,k) = P(i,k) - factor*E(i,k);  
		A(k,i) = A(i,k); 
	    }
	    for (i=1; i<=dim; ++i)
	    for (k=1; k<=dim; ++k) b[i] += A(i,k) * uPrevOnNewMesh[nodes[k]]; 
	}
    }


    // --		the stiffness matrix:  

    factor = tau*theta;

    for (i=1; i<=dim; ++i)
    for (k=1; k<=i;   ++k)  A(i,k) = P(i,k) + factor*E(i,k);


    // --	  	the source term:


    if (material->SourceTerm(t.Class())) 
    {
	FORALL(s,i) s[i] = 0.0;

	if (implicitEuler)     elem.assembleSource(t,s,Jac,pattern,time);
	else if (midPointRule) elem.assembleSource(t,s,Jac,pattern,time-0.5*tau);
	else
	{
	    FORALL(s0,i) s0[i] = 0.0;
	    
	    elem.assembleSource(t,s0,Jac,pattern,time-tau);
	    elem.assembleSource(t,s, Jac,pattern,time);
	    
	    factor = 1.0-theta;
	    FORALL(s,i) s[i] = theta*s[i] + factor*s0[i];
	}

	FORALL(s,i) b[i] += tau*s[i];
    }


    // --	  	the boundary conditions:

    if (t.onBoundary())
    {
	Bool cauchyBCs = False;
	Bool neumBCs   = False;

	for (i=1; i<=dim; ++i) 
	{
	    s[i] = 0.0;
	    for (k=1; k<=i; ++k) P(i,k) =  E(i,k) = 0.0;
	}
       
	if (implicitEuler)
	{
	    cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern,time);
	    neumBCs   = elem.assembleNeumannBCs(t,s,pattern, time);
	}
	else if (midPointRule)
	{
	    Real midTime = time - 0.5*tau;
	    cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern,midTime);
	    neumBCs   = elem.assembleNeumannBCs(t,s,pattern, midTime);
	}
	else
	{
	    FORALL(s0,i) s0[i] = 0.0;
       
	    neumBCs   = elem.assembleNeumannBCs(t,s0,pattern,time-tau);
	    neumBCs   = elem.assembleNeumannBCs(t,s,pattern, time);
       
	    cauchyBCs = elem.assembleCauchyBCs(t,E,s0,pattern,time-tau);
	    cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern, time);
       
	    factor = 1.0-theta;

	    if (cauchyBCs || neumBCs) s[i] = theta*s[i] + factor*s0[i] ;
	    if (cauchyBCs)
	    {
		for (i=1; i<=dim; ++i)
		for (k=1; k<=i;   ++k)  P(i,k) = theta*P(i,k) + factor*E(i,k);
	    }
	}

	if (cauchyBCs || neumBCs)  FORALL(b,i) b[i] += tau*s[i];
	if (cauchyBCs)
	{
	    for (i=1; i<=dim; ++i) 
	    for (k=1; k<=i; ++k)  A(i,k) += tau*P(i,k);
	}
    }


    if (Mesh()->innerBoundary) 
    {
	FORALL(s,i) s[i] = 0.0;

	if (implicitEuler)     elem.assembleInnerBCs(t,s,pattern,time);
	else if (midPointRule) elem.assembleInnerBCs(t,s,pattern,time-0.5*tau);
	else
	{
	    FORALL(s0,i) s0[i] = 0.0;
	    
	    elem.assembleInnerBCs(t,s0,pattern,time-tau);
	    elem.assembleInnerBCs(t,s, pattern,time);
	    
	    factor = 1.0-theta;
	    FORALL(s,i) s[i] = theta*s[i] + factor*s0[i];
	}

	FORALL(s,i) b[i] += tau*s[i];
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// --	assemble right-hand side for time error estimation and 
//					multiplicative defect correction 


void TransientHeatConduction:: assembleGlobalDefect(MLMatrix& AD, 
						    Vector<Num>& MD, 
						    Vector<Num>& bD)
{
    int    i, node;
    PATCH* patch;
    const int dim = element->NoOfNodes();
    Matrix<Num> AElem(dim,dim);
    Vector<Num> bElem(dim), MElem(dim);
    Vector<int> nodes(dim);
    Jacobian Jac;

    AD.reset();
    FORALL(bD,i) bD[i] = MD[i] = 0.0;
 
    Mesh()->resetElemIter();
    while ((patch = Mesh()->elemIterAll()))
    {
	patch->compJinv(Jac);
	element->initAb(AElem,bElem);
	  
	assembleDefect(*element, *patch, Jac, AElem, MElem, bElem);
  
	interface->getGlobalNodes(patch, nodes);

	AD.store(AElem, nodes);

	FORALL(nodes,i) 
	  if ((node=nodes[i]))
	  {
	      bD[node] += bElem[i];
	      MD[node] += MElem[i];
	  }
    }
}
//-------------------------------------------------------------------------


void TransientHeatConduction:: assembleDefect(const Element& elem, 
					      const PATCH& t, 
					      const Jacobian& Jac, 
					      Matrix<Num>& AD, Vector<Num>& MD,
					      Vector<Num>& bD)
{
    int   i, k;
    const int dim = elem.NoOfNodes();
    Matrix<Real> P(dim,dim);		 // 'parabolic' sub-matrix (mass-matrix)
    Vector<Real> s0(dim);

    for (i=1; i<=dim; ++i)
    for (k=1; k<=i;   ++k)  P(i,k) = 0.0;

    elem.assembleEllip(t,AD,Jac,0);

    if (lumpedMass) 
    {
	elem.assembleLumpedP(t,P,Jac,0);
	FORALL(MD,i) MD[i] = P(i,i);
    }
    else
    {
	elem.assembleP(t,P,Jac,0);

	for (i=1; i<=dim; ++i)
	for (k=1; k< i;   ++k) P(k,i) = P(i,k);

	for (i=1; i<=dim; ++i)			// lump Mass by row sum
	{
	    MD[i] = 0.0;
	    for (k=1; k<=dim; ++k) MD[i] += P(i,k);
	}
    }

    if (material->SourceTerm(t.Class())) 
    {
	FORALL(s0,i) s0[i] = 0.0;

	elem.assembleSource(t,bD,Jac,0,time);
	elem.assembleSource(t,s0,Jac,0,time-tau);

	FORALL(bD,i) bD[i] -= s0[i];
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

//#include "feplot1.h"	
//#include "triang1.h"	


void TransientProblem:: timeErrorByMidPointRule(Real* timeError, 
						Real* UNorm, Real TOL)
{
    int  i;
    
    *UNorm = compUNorm(u);

    Vector<Num> u2(u.high());
    FORALL(u,i) u2[i] = u[i];
 
    const Bool midPointRule0 = midPointRule;
    const Real theta0 = theta;
    midPointRule  = True;
    theta = 0.5;

    Problem::assembleGlobal(time);

    midPointRule = midPointRule0;
    theta = theta0;

    Real precision = sqr(rhoSpace*TOL);
    precision *= quot(sqr(*UNorm),Energy.Top());
    int step = Energy.high();

    Ab->solve(u2, 0.0, 0.0, precision, step);

    FORALL(u,i) u2[i] -= u[i];
    *timeError = compUNorm(u2);
}
//-------------------------------------------------------------------------


// --	this routine uses a multiplicative defect correction scheme 
//						 for time error estimation


void TransientProblem:: compTimeError(Real* timeError, Real* UNorm, Real TOL)
{
    if (Cmd.isSet("TError", "midPointRule"))		// for test purposes
    {
	timeErrorByMidPointRule(timeError, UNorm, TOL);
	return;
     }

    static Real accTime = 0.0;
    Timer timer, accTimer;

    int  i;
    const int dim = u.high();

    MLMatrix* AD;
    if (nComp == 1) AD = new MLSparseMatrix(sym, SpaceDim(), dim);
    else	    AD = new MLBlockMatrix (sym, SpaceDim(), dim, nComp);

    Vector<Num> u2(dim), MD(dim), bD(dim), eta(dim);


    assembleGlobalDefect(*AD,MD,bD);

    *UNorm = compUNorm(u, &MD);

    FORALL(bD,i) u2[i] = u[i] - uPrevOnNewMesh[i];

    AD->Mult(eta,u2);

    FORALL(bD,i) bD[i] = 0.5*tau * (eta[i] - bD[i]);  	  // new right hand side
    FORALL(bD,i) if (dirichletBCs->isSet(i)) bD[i] = 0.0;
    dirichletBCs->setValuesToZero();

    Ab->setRhs(bD);

    Real precision = sqr(rhoSpace*TOL);
    precision *= quot(sqr(*UNorm) , Energy.Top());

    int step = Energy.high();
    FORALL(eta,i) eta[i] = 0.0;

    Ab->solve(eta, 0.0, 0.0, precision, step);

    *timeError = compUNorm(eta, &MD);


    if (Cmd.isTrue("timeTimeError"))
    	 { cout << "\n\tTimeError Estimation: "; timer.cpu(); }
    if (Cmd.isTrue("accTimeTimeError"))
    {
	accTime += accTimer.cpu(False);
	cout << Form("\tAccumulated time:  %1.2f sec.\n", accTime); 
    }

    delete AD;

    //-------------------------  test output ------------------------

    /*
       MESH1 *mesh = (MESH1*) Mesh();
       FEPlotMESH1 Plot(mesh, X11, "DefectCorrection", 0.4);
       Plot.plotElements();
       Plot.plotBoundary();
       FORALL(bD,i) u2[i] = u[i] + eta[i];
       Plot.plotSolution(eta);
       Pause();
       
       Real eNorm = sqrt(2.0*EnergyNorm(eta));
       cout << "\n--- Time Error: eNorm=" << eNorm << "  L2Mass= " 
       << *timeError << "  E/M = " << quot(eNorm,(*timeError)) << "\n";
       
       AD.Mult(bD,eta);
       eNorm = sqrt(cdot(bD,eta));
       Real tError = sqrt(2.0*L2MassNorm(eta));
       cout << "--- Time Error: ellipNorm=" << eNorm << "  L2Mass= " << tError
       << "  Ell/M = " << quot(eNorm,tError) << "\n";
       
       AD.Mult(bD,u);
       eNorm = sqrt(2.0*EllipticNorm(u));
       cout << "--- Time Sol: eNorm=" << eNorm << "  UNorm= " << *UNorm
       << "  E/M = " << quot(eNorm,(*UNorm)) << "\n";
       */
}