[occt.git] / src / FEmTool / FEmTool_Assembly.cxx

// Created on: 1998-11-17
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


#include <FEmTool_Assembly.ixx>
#include <FEmTool_ListIteratorOfListOfVectors.hxx>
#include <FEmTool_ListOfVectors.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <math_Matrix.hxx>

//----------------------------------------------------------------------------
// Purpose - to find min index of global variables and define
//----------------------------------------------------------------------------
static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imin;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imin = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imin = Min(Imin, T->Value(nvar));
      }
    }
  return Imin;
}

//----------------------------------------------------------------------------
// Purpose - to find max index of global variables 
//----------------------------------------------------------------------------
static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imax;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imax = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imax = Max(Imax, T->Value(nvar));
      }
    }
  return Imax;
}


//=======================================================================
//function : FEmTool_Assembly
//purpose  : 
//=======================================================================
FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
				   const Handle(FEmTool_HAssemblyTable)& Table):
     myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
     B(MinIndex(Table), MaxIndex(Table))
     
{
  IsSolved = Standard_False;
  myDepTable = Dependence;
  myRefTable = Table;

  TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());

  Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;

  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T;
  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {

      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      Imin = T->Value(nvarl) + I0;
     
      for(nvar = nvarl; nvar <= nvaru; nvar++) 
	Imin = Min(Imin, T->Value(nvar) + I0);
      
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	i = T->Value(nvar) + I0;
	FirstIndexes(i) = Min(FirstIndexes(i), Imin);
      }
    }


  H = new FEmTool_ProfileMatrix(FirstIndexes);

  NullifyMatrix();
  NullifyVector();
}

void FEmTool_Assembly::NullifyMatrix() 
{
  H->Init(0.);
  IsSolved = Standard_False;
}


//=======================================================================
//function : AddMatrix
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
				 const Standard_Integer Dimension1,
				 const Standard_Integer Dimension2,
				 const math_Matrix& Mat)
{

  if(myDepTable(Dimension1, Dimension2) == 0) 
    Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");

  const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
  const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();

  Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);


  Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,

                   i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;

  for(i = nvarl; i <= nvaru; i++) {
    I = T1(i) + I0;
    ii = i0+i;
    for(j = nvarl; j <= i; j++) {
      J = T2(j) + I0;
      H->ChangeValue(I, J) += Mat(ii, j0+j);
    }
  }

  IsSolved = Standard_False;
}


//=======================================================================
//function : NullifyVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::NullifyVector() 
{
  B.Init(0.);
}

//=======================================================================
//function : AddVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddVector(const Standard_Integer Element,
				 const Standard_Integer Dimension,
				 const math_Vector& Vec) 
{
  const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
  Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
                   i0 = Vec.Lower() - nvarl;

//  Standard_Integer I, i;
  Standard_Integer i;
  
  for(i = nvarl; i <= nvaru; i++) 
    B(T(i)) += Vec(i0 + i);
  
}


//=======================================================================
//function : Solve
//purpose  : 
//=======================================================================
Standard_Boolean FEmTool_Assembly::Solve() 
{
  IsSolved = H->Decompose();

#if DEB
  if (!IsSolved) {
    cout << "Solve Echec  H = " << endl;
    H->OutM();
  }
#endif

/* ????
  while(!IsSolved && count < 5) {
    Standard_Integer i;
    for(i = 1; i <= H->RowNumber(); i++) {
      H->ChangeValue(i,i) *= 2.;
    }
    IsSolved = H->Decompose();
    count++;
  }
*/    

// Standard_Integer count = 0;
  if(!G.IsEmpty() && IsSolved) {
    // calculating H-1 Gt
    math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());

    if(GHGt.IsNull() || GHGt->RowNumber() != G.Length()) {
      TColStd_Array1OfInteger FirstIndexes(1, G.Length());
//-----------------------------------------------------------------
      TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
      Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
      Standard_Boolean Block, Zero;
      for(i = 2; i <= NbGlobVar(); i++) {
	BlockEnd = i - 1;
	if(!H->IsInProfile(i, BlockEnd)) {
	  // Maybe, begin of block
	  Block = Standard_True;
	  for(j = i + 1; j <= NbGlobVar(); j++) {
	    if(H->IsInProfile(j, BlockEnd)) {
	      Block = Standard_False;
	      break;
	    }
	  }
	  if(Block) {
	    for(j = i; j <= NbGlobVar(); j++) {
	      for(k = BlockBeg; k <= BlockEnd; k++) {
		H1(j, k) = 0; H1(k, j) = 0;
	      } 
	    }
	    BlockBeg = BlockEnd + 1;
	  }
	  else i = j;
	}
      }
      
      FEmTool_ListIteratorOfListOfVectors Iter1;
      FEmTool_ListIteratorOfListOfVectors Iter2;
      for(i = 1; i <= G.Length(); i++) {
	const FEmTool_ListOfVectors& Gi = G.Value(i);
	for(j = 1; j <= i; j++) {
	  const FEmTool_ListOfVectors& Gj = G.Value(j);
	  Zero = Standard_True;
	  for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
	    const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
	    for(k = a->Lower(); k <= a->Upper(); k++) {
	      for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
		const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
		for(l = b->Lower(); l <= b->Upper(); l++) {
		  if(H1(k, l) != 0) {
		    Zero = Standard_False;
		    break;
		  }
		}
		if(!Zero) break;
	      }
	      if(!Zero) break;
	    }
	    if(!Zero) break;
	  }
	  if(!Zero) {
	    FirstIndexes(i) = j;
	    break;
	  }
	}
      }
//-----------------------------------------------------------------------
//      for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
//	cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
      //      FirstIndexes.Init(1); // temporary GHGt is full matrix
      GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
    }

    GHGt->Init(0.);
    Standard_Integer i, j, k;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {

      const FEmTool_ListOfVectors& L = G.Value(i);
      gi.Init(0.);
      // preparing i-th line of G (or column of Gt) 
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();

	for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G 

      }
                        //                                     -1 t
      H->Solve(gi, qi); // solving H*qi = gi, qi is column of H  G


      //                               -1 t
      // Calculation of product M = G H G
      // for each i all elements of i-th column of M are calculated for k >= i
      for(k = i; k <= G.Length(); k++) {

	if(GHGt->IsInProfile(k, i)) {
	  Standard_Real m = 0.; // m = M(k,i)
	
	  const FEmTool_ListOfVectors& L = G.Value(k);
	
	  for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	    const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	    for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
	                                             // k-th line of G and i-th column of H-1 Gt 

	  }
	
	  GHGt->ChangeValue(k, i) = m;
	}
      }

    }


    IsSolved = GHGt->Decompose();
/*    count = 0;
    while(!IsSolved && count < 5) {
      for(i = 1; i <= GHGt->RowNumber(); i++) {
	GHGt->ChangeValue(i,i) *= 2.;
      }
      IsSolved = GHGt->Decompose();
      count++;
    }*/

  }

  return IsSolved;

}


//=======================================================================
//function : Solution
//purpose  : 
//=======================================================================
void FEmTool_Assembly::Solution(math_Vector& Solution) const
{
  if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");

  if(G.IsEmpty()) H->Solve(B, Solution);
  else {

    math_Vector v1(B.Lower(), B.Upper());
    H->Solve(B, v1);
    
    math_Vector l(1, G.Length()), v2(1, G.Length());
    Standard_Integer i, j;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {
      
      const FEmTool_ListOfVectors& L = G.Value(i);
      Standard_Real m = 0.;
      
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	
	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	for(j = a->Lower(); j <= a->Upper(); j++) 
	  m += v1(j) * a->Value(j); // scalar product
	// G v1
      }
      
      v2(i) = m - C.Value(i);
    }
    
    GHGt->Solve(v2, l); // Solving M*l = v2
    
    v1 = B;
    // Calculation v1 = B-Gt*l
    // v1(j) = B(j) - Gt(j,i)*l(i) = B(j) - G(i,j)*l(i)
    // 
      for(i = 1; i <= G.Length(); i++) {
	
	const FEmTool_ListOfVectors& L = G.Value(i);
	
	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	  
	  const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	  for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
	  
	}
	
      }
      
    H->Solve(v1, Solution);
  }
  
}

Standard_Integer FEmTool_Assembly::NbGlobVar() const
{

  return B.Length();

}

void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
{
  AssTable = myRefTable;
}


void FEmTool_Assembly::ResetConstraint() 
{
  G.Clear();
  C.Clear();
}

void FEmTool_Assembly::NullifyConstraint() 
{
  FEmTool_ListIteratorOfListOfVectors Iter;
  Standard_Integer i;

  for(i = 1; i <= G.Length(); i++) {
  
    C.SetValue(i, 0.);

    for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next()) 
      Iter.Value()->Init(0.);

  }

}


//=======================================================================
//function : AddConstraint
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
				     const Standard_Integer Element,
				     const Standard_Integer Dimension,
				     const math_Vector& LinearForm,
				     const Standard_Real Value) 
{
  while(G.Length() < IndexofConstraint) {
    // Add new lines in G
    FEmTool_ListOfVectors L;
    G.Append(L);
    C.Append(0.);
  }

  FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
  
  Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
  Standard_Integer i, Imax = 0, Imin = NbGlobVar();

  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
    Imin = Min(Imin, Indexes->Value(i));
    Imax = Max(Imax, Indexes->Value(i));
  }

  Handle(TColStd_HArray1OfReal) Coeff;

  if(L.IsEmpty()) {
    Coeff = new TColStd_HArray1OfReal(Imin,Imax);
    Coeff->Init(0.);
    L.Append(Coeff);
  }
  else {
    FEmTool_ListIteratorOfListOfVectors Iter(L);
    Standard_Integer i;
    Standard_Real  s1 = 0, s2 = 0;
    Handle(TColStd_HArray1OfReal) Aux1, Aux2;
    for(i=1; Iter.More(); Iter.Next(), i++) {
      if(Imin >= Iter.Value()->Lower()) {
	s1 = i;
	Aux1 = Iter.Value();
	if(Imax <= Iter.Value()->Upper()) {
	  s2 = s1;
	  Coeff = Iter.Value();
	  break;
	}
      }

      if(Imax <= Iter.Value()->Upper()) {
	s2 = i;
	Aux2 = Iter.Value();
      }
    }

    if(s1 != s2) {
      if(s1 == 0) {
	if(Imax < Aux2->Lower()) {
	  // inserting before first segment
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Prepend(Coeff);
	}
	else {
	  // merge new and first segment
	  Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	  for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	  for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	  L.First() = Coeff;
	}
      }
      else if(s2 == 0) {
	if(Imin > Aux1->Upper()) {
	  // append new
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Append(Coeff);
	}
	else {
	  // merge new and last segment
	  Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	  for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	  for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	  L.Last() = Coeff;
	}
      }
      else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
	// merge s1 and new
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
      }
      else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) { 
	// merge new and first segment
	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s2; Iter.Next(), i++);
	Iter.Value() = Coeff;
     }
      else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
	// inserting new between s1 and s2
	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	Coeff->Init(0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	L.InsertAfter(Coeff,Iter);
      }
      else {
	// merge s1, new, s2 and remove s2 
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
	Iter.Next();
	L.Remove(Iter);
      }
    }
  }

  // adding 
  Standard_Integer j = LinearForm.Lower();
  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
    Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
  }
  
  C.ChangeValue(IndexofConstraint) += Value;


}
Commit	Line	Data
b311480e	1	// Created on: 1998-11-17
	2	// Created by: Igor FEOKTISTOV
	3	// Copyright (c) 1998-1999 Matra Datavision
	4	// Copyright (c) 1999-2012 OPEN CASCADE SAS
	5	//
	6	// The content of this file is subject to the Open CASCADE Technology Public
	7	// License Version 6.5 (the "License"). You may not use the content of this file
	8	// except in compliance with the License. Please obtain a copy of the License
	9	// at http://www.opencascade.org and read it completely before using this file.
	10	//
	11	// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
	12	// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
	13	//
	14	// The Original Code and all software distributed under the License is
	15	// distributed on an "AS IS" basis, without warranty of any kind, and the
	16	// Initial Developer hereby disclaims all such warranties, including without
	17	// limitation, any warranties of merchantability, fitness for a particular
	18	// purpose or non-infringement. Please see the License for the specific terms
	19	// and conditions governing the rights and limitations under the License.
	20
7fd59977	21
	22
	23	#include <FEmTool_Assembly.ixx>
	24	#include <FEmTool_ListIteratorOfListOfVectors.hxx>
	25	#include <FEmTool_ListOfVectors.hxx>
	26	#include <TColStd_HArray1OfReal.hxx>
	27	#include <math_Matrix.hxx>
	28
	29	//----------------------------------------------------------------------------
	30	// Purpose - to find min index of global variables and define
	31	//----------------------------------------------------------------------------
	32	static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
	33	{
	34	Standard_Integer dim, el, nvar, Imin;
	35	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	36	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	37
	38	Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
	39
	40	nvarl = T->Lower();
	41
	42	Imin = T->Value(nvarl);
	43
	44	for(dim = diml; dim <= dimu; dim++)
	45	for(el = ell; el <= elu; el++) {
	46	T = Table->Value(dim, el);
	47	nvarl = T->Lower(); nvaru = T->Upper();
	48	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	49	Imin = Min(Imin, T->Value(nvar));
	50	}
	51	}
	52	return Imin;
	53	}
	54
	55	//----------------------------------------------------------------------------
	56	// Purpose - to find max index of global variables
	57	//----------------------------------------------------------------------------
	58	static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
	59	{
	60	Standard_Integer dim, el, nvar, Imax;
	61	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	62	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	63
	64	Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
	65
	66	nvarl = T->Lower();
	67
	68	Imax = T->Value(nvarl);
	69
	70	for(dim = diml; dim <= dimu; dim++)
	71	for(el = ell; el <= elu; el++) {
	72	T = Table->Value(dim, el);
	73	nvarl = T->Lower(); nvaru = T->Upper();
	74	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	75	Imax = Max(Imax, T->Value(nvar));
	76	}
	77	}
	78	return Imax;
	79	}
	80
	81
	82
	83	//=======================================================================
	84	//function : FEmTool_Assembly
85	//purpose :
86	//=======================================================================
87	FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
88	const Handle(FEmTool_HAssemblyTable)& Table):
89	myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
90	B(MinIndex(Table), MaxIndex(Table))
91
92	{
93	IsSolved = Standard_False;
94	myDepTable = Dependence;
95	myRefTable = Table;
96
97	TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());
98
7fd59977	99	Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;
6e6cd5d9	100
7fd59977	101	Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
	102	ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
	103
	104	Handle(TColStd_HArray1OfInteger) T;
	105	for(dim = diml; dim <= dimu; dim++)
	106	for(el = ell; el <= elu; el++) {
	107
	108	T = Table->Value(dim, el);
	109	nvarl = T->Lower(); nvaru = T->Upper();
	110	Imin = T->Value(nvarl) + I0;
	111
	112	for(nvar = nvarl; nvar <= nvaru; nvar++)
	113	Imin = Min(Imin, T->Value(nvar) + I0);
	114
	115	for(nvar = nvarl; nvar <= nvaru; nvar++) {
	116	i = T->Value(nvar) + I0;
	117	FirstIndexes(i) = Min(FirstIndexes(i), Imin);
	118	}
	119	}
	120
	121
	122	H = new FEmTool_ProfileMatrix(FirstIndexes);
	123
	124	NullifyMatrix();
	125	NullifyVector();
	126	}
	127
	128	void FEmTool_Assembly::NullifyMatrix()
	129	{
	130	H->Init(0.);
	131	IsSolved = Standard_False;
	132	}
	133
	134
	135	//=======================================================================
	136	//function : AddMatrix
	137	//purpose :
	138	//=======================================================================
	139	void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
	140	const Standard_Integer Dimension1,
	141	const Standard_Integer Dimension2,
	142	const math_Matrix& Mat)
	143	{
	144
	145	if(myDepTable(Dimension1, Dimension2) == 0)
	146	Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");
	147
	148	const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
	149	const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();
	150
	151	Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);
	152
6e6cd5d9	153
7fd59977	154	Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,
6e6cd5d9	155
7fd59977	156	i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;
	157
	158	for(i = nvarl; i <= nvaru; i++) {
	159	I = T1(i) + I0;
	160	ii = i0+i;
	161	for(j = nvarl; j <= i; j++) {
	162	J = T2(j) + I0;
	163	H->ChangeValue(I, J) += Mat(ii, j0+j);
	164	}
	165	}
	166
	167	IsSolved = Standard_False;
	168	}
	169
	170
	171	//=======================================================================
	172	//function : NullifyVector
	173	//purpose :
	174	//=======================================================================
	175	void FEmTool_Assembly::NullifyVector()
	176	{
	177	B.Init(0.);
	178	}
	179
	180	//=======================================================================
	181	//function : AddVector
	182	//purpose :
	183	//=======================================================================
	184	void FEmTool_Assembly::AddVector(const Standard_Integer Element,
	185	const Standard_Integer Dimension,
	186	const math_Vector& Vec)
	187	{
	188	const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
	189	Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
	190	i0 = Vec.Lower() - nvarl;
	191
	192	// Standard_Integer I, i;
	193	Standard_Integer i;
	194
	195	for(i = nvarl; i <= nvaru; i++)
	196	B(T(i)) += Vec(i0 + i);
	197
	198	}
	199
	200
	201	//=======================================================================
	202	//function : Solve
	203	//purpose :
	204	//=======================================================================
	205	Standard_Boolean FEmTool_Assembly::Solve()
	206	{
	207	IsSolved = H->Decompose();
	208
	209	#if DEB
	210	if (!IsSolved) {
	211	cout << "Solve Echec H = " << endl;
	212	H->OutM();
	213	}
	214	#endif
	215
	216	/* ????
	217	while(!IsSolved && count < 5) {
	218	Standard_Integer i;
	219	for(i = 1; i <= H->RowNumber(); i++) {
220	H->ChangeValue(i,i) *= 2.;
221	}
222	IsSolved = H->Decompose();
223	count++;
224	}
225	*/
226
227	// Standard_Integer count = 0;
228	if(!G.IsEmpty() && IsSolved) {
229	// calculating H-1 Gt
230	math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());
231
232	if(GHGt.IsNull() \|\| GHGt->RowNumber() != G.Length()) {
233	TColStd_Array1OfInteger FirstIndexes(1, G.Length());
234	//-----------------------------------------------------------------
235	TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
236	Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
237	Standard_Boolean Block, Zero;
238	for(i = 2; i <= NbGlobVar(); i++) {
239	BlockEnd = i - 1;
240	if(!H->IsInProfile(i, BlockEnd)) {
241	// Maybe, begin of block
242	Block = Standard_True;
243	for(j = i + 1; j <= NbGlobVar(); j++) {
244	if(H->IsInProfile(j, BlockEnd)) {
245	Block = Standard_False;
246	break;
247	}
248	}
249	if(Block) {
250	for(j = i; j <= NbGlobVar(); j++) {
251	for(k = BlockBeg; k <= BlockEnd; k++) {
252	H1(j, k) = 0; H1(k, j) = 0;
253	}
254	}
255	BlockBeg = BlockEnd + 1;
256	}
257	else i = j;
258	}
259	}
260
261	FEmTool_ListIteratorOfListOfVectors Iter1;
262	FEmTool_ListIteratorOfListOfVectors Iter2;
263	for(i = 1; i <= G.Length(); i++) {
264	const FEmTool_ListOfVectors& Gi = G.Value(i);
265	for(j = 1; j <= i; j++) {
266	const FEmTool_ListOfVectors& Gj = G.Value(j);
267	Zero = Standard_True;
268	for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
269	const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
270	for(k = a->Lower(); k <= a->Upper(); k++) {
271	for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
272	const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
273	for(l = b->Lower(); l <= b->Upper(); l++) {
274	if(H1(k, l) != 0) {
275	Zero = Standard_False;
276	break;
277	}
278	}
279	if(!Zero) break;
280	}
281	if(!Zero) break;
282	}
283	if(!Zero) break;
284	}
285	if(!Zero) {
286	FirstIndexes(i) = j;
287	break;
288	}
289	}
290	}
291	//-----------------------------------------------------------------------
292	// for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
293	// cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
294	// FirstIndexes.Init(1); // temporary GHGt is full matrix
295	GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
296	}
297
298	GHGt->Init(0.);
299	Standard_Integer i, j, k;
300	FEmTool_ListIteratorOfListOfVectors Iter;
301
302	for(i = 1; i <= G.Length(); i++) {
303
304	const FEmTool_ListOfVectors& L = G.Value(i);
305	gi.Init(0.);
306	// preparing i-th line of G (or column of Gt)
307	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
308
309	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
310
311	for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G
312
313	}
314	// -1 t
315	H->Solve(gi, qi); // solving H*qi = gi, qi is column of H G
316
317
318	// -1 t
319	// Calculation of product M = G H G
320	// for each i all elements of i-th column of M are calculated for k >= i
321	for(k = i; k <= G.Length(); k++) {
322
323	if(GHGt->IsInProfile(k, i)) {
324	Standard_Real m = 0.; // m = M(k,i)
325
326	const FEmTool_ListOfVectors& L = G.Value(k);
327
328	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
329
330	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
331	for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
332	// k-th line of G and i-th column of H-1 Gt
333
334	}
335
336	GHGt->ChangeValue(k, i) = m;
337	}
338	}
339
340	}
341
342
343	IsSolved = GHGt->Decompose();
344	/* count = 0;
345	while(!IsSolved && count < 5) {
346	for(i = 1; i <= GHGt->RowNumber(); i++) {
347	GHGt->ChangeValue(i,i) *= 2.;
348	}
349	IsSolved = GHGt->Decompose();
350	count++;
351	}*/
352
353	}
354
355	return IsSolved;
356
357	}
358
359
360	//=======================================================================
361	//function : Solution
362	//purpose :
363	//=======================================================================
364	void FEmTool_Assembly::Solution(math_Vector& Solution) const
365	{
366	if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");
367
368	if(G.IsEmpty()) H->Solve(B, Solution);
369	else {
370
371	math_Vector v1(B.Lower(), B.Upper());
372	H->Solve(B, v1);
373
374	math_Vector l(1, G.Length()), v2(1, G.Length());
375	Standard_Integer i, j;
376	FEmTool_ListIteratorOfListOfVectors Iter;
377
378	for(i = 1; i <= G.Length(); i++) {
379
380	const FEmTool_ListOfVectors& L = G.Value(i);
381	Standard_Real m = 0.;
382
383	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
384
385	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
386	for(j = a->Lower(); j <= a->Upper(); j++)
387	m += v1(j) * a->Value(j); // scalar product
388	// G v1
389	}
390
391	v2(i) = m - C.Value(i);
392	}
393
394	GHGt->Solve(v2, l); // Solving M*l = v2
395
396	v1 = B;
397	// Calculation v1 = B-Gt*l
398	// v1(j) = B(j) - Gt(j,i)l(i) = B(j) - G(i,j)l(i)
399	//
400	for(i = 1; i <= G.Length(); i++) {
401
402	const FEmTool_ListOfVectors& L = G.Value(i);
403
404	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
405
406	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
407	for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
408
409	}
410
411	}
412
413	H->Solve(v1, Solution);
414	}
415
416	}
417
418	Standard_Integer FEmTool_Assembly::NbGlobVar() const
419	{
420
421	return B.Length();
422
423	}
424
425	void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
426	{
427	AssTable = myRefTable;
428	}
429
430
431	void FEmTool_Assembly::ResetConstraint()
432	{
433	G.Clear();
434	C.Clear();
435	}
436
437	void FEmTool_Assembly::NullifyConstraint()
438	{
439	FEmTool_ListIteratorOfListOfVectors Iter;
440	Standard_Integer i;
441
442	for(i = 1; i <= G.Length(); i++) {
443
444	C.SetValue(i, 0.);
445
446	for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next())
447	Iter.Value()->Init(0.);
448
449	}
450
451	}
452
453
454	//=======================================================================
455	//function : AddConstraint
456	//purpose :
457	//=======================================================================
458	void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
459	const Standard_Integer Element,
460	const Standard_Integer Dimension,
461	const math_Vector& LinearForm,
462	const Standard_Real Value)
463	{
464	while(G.Length() < IndexofConstraint) {
465	// Add new lines in G
466	FEmTool_ListOfVectors L;
467	G.Append(L);
468	C.Append(0.);
469	}
470
471	FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
472
473	Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
474	Standard_Integer i, Imax = 0, Imin = NbGlobVar();
475
476	for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
477	Imin = Min(Imin, Indexes->Value(i));
478	Imax = Max(Imax, Indexes->Value(i));
479	}
480
481	Handle(TColStd_HArray1OfReal) Coeff;
482
483	if(L.IsEmpty()) {
484	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
485	Coeff->Init(0.);
486	L.Append(Coeff);
487	}
488	else {
489	FEmTool_ListIteratorOfListOfVectors Iter(L);
490	Standard_Integer i;
491	Standard_Real s1 = 0, s2 = 0;
492	Handle(TColStd_HArray1OfReal) Aux1, Aux2;
493	for(i=1; Iter.More(); Iter.Next(), i++) {
494	if(Imin >= Iter.Value()->Lower()) {
495	s1 = i;
496	Aux1 = Iter.Value();
497	if(Imax <= Iter.Value()->Upper()) {
498	s2 = s1;
499	Coeff = Iter.Value();
500	break;
501	}
502	}
503
504	if(Imax <= Iter.Value()->Upper()) {
505	s2 = i;
506	Aux2 = Iter.Value();
507	}
508	}
509
510	if(s1 != s2) {
511	if(s1 == 0) {
512	if(Imax < Aux2->Lower()) {
513	// inserting before first segment
514	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
515	Coeff->Init(0.);
516	L.Prepend(Coeff);
517	}
518	else {
519	// merge new and first segment
520	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
521	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
522	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
523	L.First() = Coeff;
524	}
525	}
526	else if(s2 == 0) {
527	if(Imin > Aux1->Upper()) {
528	// append new
529	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
530	Coeff->Init(0.);
531	L.Append(Coeff);
532	}
533	else {
534	// merge new and last segment
535	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
536	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
537	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
538	L.Last() = Coeff;
539	}
540	}
541	else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
542	// merge s1 and new
543	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
544	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
545	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
546	Iter.Initialize(L);
547	for(i = 1; i < s1; Iter.Next(), i++);
548	Iter.Value() = Coeff;
549	}
550	else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) {
551	// merge new and first segment
552	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
553	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
554	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
555	Iter.Initialize(L);
556	for(i = 1; i < s2; Iter.Next(), i++);
557	Iter.Value() = Coeff;
558	}
559	else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
560	// inserting new between s1 and s2
561	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
562	Coeff->Init(0.);
563	Iter.Initialize(L);
564	for(i = 1; i < s1; Iter.Next(), i++);
565	L.InsertAfter(Coeff,Iter);
566	}
567	else {
568	// merge s1, new, s2 and remove s2
569	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
570	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
571	for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
572	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
573	Iter.Initialize(L);
574	for(i = 1; i < s1; Iter.Next(), i++);
575	Iter.Value() = Coeff;
576	Iter.Next();
577	L.Remove(Iter);
578	}
579	}
580	}
581
582	// adding
583	Standard_Integer j = LinearForm.Lower();
584	for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
585	Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
586	}
587
588	C.ChangeValue(IndexofConstraint) += Value;
589
590
591
592	}
593