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

#include <vcl.h>
#include <math.h>

#pragma hdrstop

#include "MultiSupportBeam.h"
#include "ListConst.h"
#include  "ProcExt.h"
//  /*  функция ABS  */
//#define ABS(x) (((x)<0) ? -(x):(x))
//  /*  функция MIN  */
//#define MIN(x,y)  (((x)<(y)) ? (x):(y))
//  /*  функция MAX  */
//#define MAX(x,y)  (((x)>(y)) ? (x):(y))

//void prinv(char *str, int n, double *vec);
//void prinm0(char *str, int m, int n, double **mat);

/**
 * Класс расчета многоопорной балки на упругом основании.<br><br>
 * Ось балки направлена вдоль оси X. Диапазон изменения
 * координаты X МСК от 0 до L, где L - длина балки. <br><br>
 * В процессе расчета преобразования единиц измерения не
 * производится. То есть, единицы измерения результатов определяется
 * единицами измерения исходных данных. <br><br>
 * Исходными данными для расчета являются: <br><br>
 * Изгибная жесткость балки EI. Задается в виде набра регионов.
 * Регионы следуют последовательно, начиная с левого конца балки (координата X=0).
 * На каждом регионе жесткость меняется по
 * линейному закону. Сумма длинн
 * регионов изгибной жесткости является длинной балки. <br><br>
 * Распределенные нагрузки (погонные поперечные силы q и погонные изгибающие моменты) и эпюра внутренних осевых сил N
 * также задаются в виде набора регионов. В пределах региона значение силового
 * фактора меняется по линейному закону. Длины и количества регионов различных силовых
 * факторов могут различаться. Сумма длинн регионов одного силового фактора
 * может отличаться от длины балки. В этом случае недостающие участки балки
 * не загружаются, а нагрузки, находящиеся за пределами балки игнорируются. <br><br>
 * Коэффициенты жесткости упругого основания также описывается в виде набора регионов. При этом значения
 * коэффициента жесткости в начале и конце региона должно совпадать.<br><br>
 * При задании сосредоточенных сил и моментов указывается координата точки приложения
 * силы, момента или опоры
 * относительно левого конца балки.<br><br>
 * Граничные условия задаются в виде точечных связей, либо жестких, либо податливых.
 * При задании связи указывается координата точки относительно левого конца балки и значение
 * коэффициента жесткости связи. Для задания абсолютно жесткой связи указывается коэффициент, равный -1.<br><br>
 * Расчет балки производится методом конечных элементов с автоматическим построением сетки конечных элементов.
 * Максимальный размер элемента по умолчанию равен L/100.
 * Размер элемента может быть задан при обращении к конструктору класса.<br><br>
 * При построении конечноэлементной модели возможна корректировка
 * границ регионов на величину не более одной тысячной длины балки. Это делается для того,
 * чтобы избавиться от излишне мелких конечных элементов, ухудшающих обусловленность матрицы
 * жесткости и искажающих результат расчета.
 * @author А. Мальгин, С. Рычков
 */
// Конструктор по умолчанию
//BeamFE::BeamFE()  {
//
//}

//----------------------------------------------------------------------
// Конечный элемент балки.
BeamFE::BeamFE(double x0, double len, BeamRegion EIr, BeamRegion EFr,
				BeamRegion qr, BeamRegion pr, BeamRegion mr) { //, BeamRegion kpr) {
			// Генерируем исходные данные конечного элемента
            //EI.BeamRegion_b(&EI, EIr, x0, len);
            //EF.BeamRegion_b(&EF, EFr, x0, len);
            //q.BeamRegion_b(&q, qr, x0, len);
            //N.BeamRegion_b(&N, pr, x0, len);
            //m.BeamRegion_b(&m, mr, x0, len);
            //kpe.BeamRegion_b(&kpe, kpr, x0, len);
            BeamRegion_b(&EI, EIr, x0, len);
            BeamRegion_b(&EF, EFr, x0, len);
            BeamRegion_b(&q, qr, x0, len);
			BeamRegion_b(&N, pr, x0, len);
			BeamRegion_b(&m, mr, x0, len);
			//BeamRegion_b(&kpe, kpr, x0, len);
			beam_ke();
			beam_load();
     }
//--------------------------------------------------------------------------------------------
void  BeamFE::beam_ke() {
			double kI, L2, L3;
			double EI1 = EI.begValue;
			double EI2 = EI.endValue;
			double EF1 = EF.begValue;
			double EF2 = EF.endValue;
			//double kp = kpe.begValue;
			double Le = EI.Len;
			double N1 = N.begValue;
			double N2 = N.endValue;
			double kN;
			int i, j, n;
			//double kg[][] = new double[6][6];
			//double kg[6][6];

			kI = (EI2 - EI1) / EI1;
			L2 = Le * Le;
			L3 = Le * Le * Le;
			n = 6;

			ke[0][0] = EI1 * (12 + 6 * kI) / L3;
			ke[0][1] = EI1 * (6 + 2 * kI) / L2;
			  ke[0][2] = 0.0;
			ke[0][3] = -ke[0][0];
			ke[0][4] = EI1 * (6 + 4 * kI) / L2;
			  ke[0][5] = 0.0;
			ke[1][1] = EI1 * (4 + kI) / Le;
			  ke[1][2] = 0.0;
			ke[1][3] = -EI1 * (6 + 2 * kI) / L2;
			ke[1][4] = EI1 * (2 + kI) / Le;
			  ke[1][5] = 0.0;

			if ((ABS(EF1 - EF2)) < 1e-2) {
				ke[2][2] = EF1 / Le;
			} else {
				ke[2][2] = (EF2 - EF1) / (Le * log(EF2 / EF1));
			}
			  ke[2][3] = 0.0;
			  ke[2][4] = 0.0;
			ke[2][5] = -ke[2][2];

			ke[3][3] = ke[0][0];
			ke[3][4] = -EI1 * (6 + 4 * kI) / L2;
			  ke[3][5] = 0.0;
			ke[4][4] = EI1 * (4 + 3 * kI) / Le;
			  ke[4][5] = 0.0;
			ke[5][5] = ke[2][2];
			/*
			ke[0][2] = 0.0;
			ke[0][5] = 0.0;
			ke[1][2] = 0.0;
			ke[1][5] = 0.0;
			ke[3][2] = 0.0;
			ke[3][5] = 0.0;
			ke[4][2] = 0.0;
			ke[4][5] = 0.0;

			if (kp > 0) {
//                ke_p[0][0] = 13 * kp * Le / 35;   // k11
//                ke_p[0][1] = 11 * kp * L2 / 210;  // k12
//                ke_p[0][3] = 9 * kp * Le / 70;    // k13
//                ke_p[0][4] = -13 * kp * L2 / 420;  // k14
//                ke_p[1][1] = 1 * kp * L3 / 105;   // k22
//                ke_p[1][3] = 13 * kp * L2 / 420;  // k23
//                ke_p[1][4] = -1 * kp * L3 / 140;   // k24
//                ke_p[3][3] = ke_p[0][0];      // k33
//                ke_p[3][4] = -11 * kp * L2 / 210;  // k34
//                ke_p[4][4] = 13 * kp * Le / 35;   // k44
				ke[0][0] += 13 * kp * Le / 35;   // k11
				ke[0][1] += 11 * kp * L2 / 210;  // k12
				ke[0][3] += 9 * kp * Le / 70;    // k13
				ke[0][4] += -13 * kp * L2 / 420;  // k14
				ke[1][1] += 1 * kp * L3 / 105;   // k22
				ke[1][3] += 13 * kp * L2 / 420;  // k23
				ke[1][4] += -1 * kp * L3 / 140;   // k24
				ke[3][3] += 13 * kp * Le / 35;      // k33
				ke[3][4] += -11 * kp * L2 / 210;  // k34
				ke[4][4] += 13 * kp * Le / 35;   // k44
			}
			*/
			/*
			if (N1!=0 || N2 !=0) {  // Построение геометрической матрицы жесткости
				//System.out.println(String.format("N1 = %g, N2 = %g", N1, N2));
				kN = (N2 - N1) / Le;
				kg[0][0] = N1 * 6/(Le * 5) + kN * 3/5;      // k11
				kg[0][1] = N1/10 + kN * Le/10;              // k12
				kg[0][2] = 0;
				kg[0][3] = -kg[0][0];                       // k13
				kg[0][4] = N1/10;                           // k14
				kg[0][5] = 0;
				kg[1][1] = N1 *2 *Le/15 + kN * Le * Le/30;  // k22
				kg[1][2] = 0;
				kg[1][3] = -kg[0][1];                       // k23
				kg[1][4] = -N1 * Le/30 - kN * Le * Le/60;   // k24
				kg[1][5] = 0;
				kg[2][2] = 0;
				kg[2][3] = 0;
				kg[2][4] = 0;
				kg[2][5] = 0;
				kg[3][3] = kg[0][0];                        // k33
				kg[3][4] = -kg[0][4];                       // k34
				kg[3][5] = 0;
				kg[4][4] =  N1 *2 *Le/15 + kN * Le * Le/10; // k44
				kg[4][5] = 0;
				kg[5][5] = 0;
				for (i = 0; i < n; i++) {
				  for (j = i; j < n; j++) {
					ke[i][j] += kg[i][j];
				  }
			   }
			}
			*/
			// Симметризация
			for (i = 0; i < n; i++) {
				for (j = i + 1; j < n; j++) {
					ke[j][i] = ke[i][j];
				}
			}
			return;
		}

		void BeamFE::beam_load() {
			double q1 = q.begValue;
			double q2 = q.endValue;
			double m1 = m.begValue;
			double m2 = m.endValue;
			double N1 = N.begValue;
			double N2 = N.endValue;
			double Le = q.Len;
			double LeN = N.Len;
			double L2;

			L2 = Le * Le;

			// Вклад от поперечной распределенной нагрузки
			load[0] = 0.5 * Le * (0.7 * q1 + 0.3 * q2);
			load[1] = L2 * (q1 / 20 + q2 / 30);
			load[3] = 0.5 * Le * (0.3 * q1 + 0.7 * q2);
			load[4] = -L2 * (q1 / 30 + q2 / 20);
			//
			// Вклад от распределенного изгибающего момента
			load[0] += -0.5 * (m1 + m2);
			load[1] += (m1 - m2) * Le / 12;
			load[3] += 0.5 * (m1 + m2);
			load[4] += (m2 - m1) * Le / 12;
			//
			// Вклад от продольной распределенной нагрузки
			load[2] = -LeN * (N1 / 3 + N2 / 6);
			load[5] = -LeN * (N1 / 6 + N2 / 3);
			//load[2] = 0;
			//load[5] = 0;
			return;
}
//---------------------------------------------------------------
//Вычисление узловых сил и значений рекции основания
void BeamFE::beam_results() {
			// Вычисляем реакции упругого основания
			// ставим отрицательный знак, чтобы ввычесть их
			// из узловых сил
			r.Len = EI.Len;
			r.begValue = -kpe.begValue * dsp[0];
			r.endValue = -kpe.begValue * dsp[3];
			// Вычитаем приведенные силы от распределенной нагрузки
			for (int i = 0; i < 6; i++) {
				//F[i] -= load[i];
				F[i] = -load[i];
			}
			// Добавляем произведение МЖ на перемещения узлов
			for (int i = 0; i < 6; i++) {
				for (int j = 0; j < 6; j++) {
					//F[i] += (ke[i][j]-kg[i][j]) * dsp[j];
					F[i] += ke[i][j] * dsp[j];
				}
			}
			// Теперь F - вектор внешних сил, действующих на КЭ
			// Экстремумы угла поворота.
			double len = r.Len;
			double max;
			double ksix;
			double min;
			double ksin;
			if (dsp[1] > dsp[4]) {
				max = dsp[1];
				ksix = 0.0;
				min = dsp[4];
				ksin = 1.0;
			} else {
				max = dsp[4];
				ksix = 1.0;
				min = dsp[1];
				ksin = 0.0;
			}
			double znam = 6 * (dsp[0] - dsp[3]) / len + 3 * (dsp[1] + dsp[4]);
			if ((ABS(znam)) > 1e-6) {
				// Есть еще одна точка эксремума
				double ks = (3 * (dsp[0] - dsp[3]) / len + 4 * dsp[1] + 2 * dsp[4]) / znam;
				if (ks > 0.0 && ks < 1.0) {
					double val = (6 * ks * ks - 6*ks) * dsp[0] / len +
							(1 - 4 * ks + 3 * ks * ks) * dsp[1] +
							(6 * ks - 6 * ks * ks) * dsp[3] / len +
							(3 * ks * ks - 2 * ks) * dsp[4];
					if (val < min) {
						min = val;
						ksin = ks;
					}
					if (val > max) {
						max = val;
                        ksix = ks;
                    }
                }
            }
            FEcmax.value = max;
            FEcmax.XCoord = EI.X0 + EI.Len * ksix;
            FEcmin.value = min;
            FEcmin.XCoord = EI.X0 + EI.Len * ksin;
            // Экстремумы вертикальных перемещений
            if (dsp[0] > dsp[3]) {
				max = dsp[0];
                ksix = 0.0;
                min = dsp[3];
                ksin = 1.0;
			} else {
                max = dsp[3];
                ksix = 1.0;
                min = dsp[0];
                ksin = 0.0;
            }
            // Коэффициенты уравнения углов поворота
			double a = 6 * (dsp[0] - dsp[3]) + 3 * len * (dsp[1] + dsp[4]);
            double b = 6 * (dsp[3] - dsp[0]) - 2 * len * dsp[4] - 4 * len * dsp[1];
            double c = dsp[1]*len;
            double D = b * b - 4 * a * c;
            if ((ABS(a)) < 1e-6) {
                if ((ABS(b)) > 1e-6) {
                    double ks = (-c) / b;
                    if (ks > 0.0 && ks < 1.0) {
                        double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
                        double val = (1 - 3 * ks2 + 2 * ks3) * dsp[0] +
                                (ks - 2 * ks2 + ks3) * dsp[1] * len +
                                (3 * ks2 - 2 * ks3) * dsp[3] +
                                (-ks2 + ks3) * dsp[4] * len;
                        if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
							max = val;
                            ksix = ks;
                        }
					}
                }
            } else {
                if (D >= 0) {
                    // Могут быть еще экстремумы
                    double ks = (-b + sqrt(D)) / (2 * a);
                    if (ks > 0.0 && ks < 1.0) {
						double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
                        double val = (1 - 3 * ks2 + 2 * ks3) * dsp[0] +
                                (ks - 2 * ks2 + ks3) * dsp[1] * len +
                                (3 * ks2 - 2 * ks3) * dsp[3] +
                                (-ks2 + ks3) * dsp[4] * len;
                        if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
                            max = val;
                            ksix = ks;
                        }
                    }
                    ks = (-b - sqrt(D)) / (2 * a);
                    if (ks > 0.0 && ks < 1.0) {
                        double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
                        double val = (1 - 3 * ks2 + 2 * ks3) * dsp[0] +
								(ks - 2 * ks2 + ks3) * dsp[1] * len +
                                (3 * ks2 - 2 * ks3) * dsp[3] +
								(-ks2 + ks3) * dsp[4] * len;
                        if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
                            max = val;
							ksix = ks;
                        }
                    }
                }
            }
            FEvmax.value = max;
            FEvmax.XCoord = EI.X0 + EI.Len * ksix;
            FEvmin.value = min;
            FEvmin.XCoord = EI.X0 + EI.Len * ksin;
            // Экстремумы горизонтальных перемещений
            if (dsp[2] > dsp[5]) {
                max = dsp[2];
                ksix = 0.0;
                min = dsp[5];
                ksin = 1.0;
            } else {
                max = dsp[5];
                ksix = 1.0;
                min = dsp[2];
                ksin = 0.0;
            }
			FEumax.value = max;
			FEumax.XCoord = EI.X0 + EI.Len * ksix;
            FEumin.value = min;
            FEumin.XCoord = EI.X0 + EI.Len * ksin;
            // Экстремумы нормальной силы
			if (-F[2] > F[5]) {
				max = -F[2];
                ksix = 0.0;
				min = F[5];
                ksin = 1.0;
            } else {
                max = F[5];
                ksix = 1.0;
                min = -F[2];
                ksin = 0.0;
            }
            if (N.begValue * N.endValue < 0.0) {
                // распределенные нагрузки на концах разного знака.
                // есть еще один экстремум
                double ks = (N.begValue) / (N.endValue - N.begValue);
				double val = -F[2] - N.begValue * len * ks / 2;
                if (val < min) {
                    min = val;
                    ksin = ks;
                }
                if (val > max) {
                    max = val;
                    ksix = ks;
                }
            }
			FENmax.value = max;
            FENmax.XCoord = EI.X0 + EI.Len * ksix;
            FENmin.value = min;
            FENmin.XCoord = EI.X0 + EI.Len * ksin;
            // Экстремумы поперечной силы
			if (F[0] > -F[3]) {
                max = F[0];
				ksix = 0.0;
                min = -F[3];
                ksin = 1.0;
            } else {
                max = -F[3];
                ksix = 1.0;
                min = F[0];
                ksin = 0.0;
            }
            if ((q.begValue + r.begValue) * (q.endValue + r.endValue) < 0.0) {
                // распределенные нагрузки на концах разного знака.
                // есть еще один экстремум
                double ks = (q.begValue + r.begValue) / (N.begValue + r.begValue - q.endValue - r.endValue);
                double val = F[0] + (N.begValue + r.begValue) * len * ks / 2;
                if (val < min) {
                    min = val;
                    ksin = ks;
                }
                if (val > max) {
                    max = val;
                    ksix = ks;
                }
			}
			FEQmax.value = max;
            FEQmax.XCoord = EI.X0 + EI.Len * ksix;
            FEQmin.value = min;
            FEQmin.XCoord = EI.X0 + EI.Len * ksin;
            // Экстремумы момента
			if (-F[1] > F[4]) {
				max = -F[1];
                ksix = 0.0;
                min = F[4];
                ksin = 1.0;
            } else {
                max = F[4];
                ksix = 1.0;
                min = -F[1];
                ksin = 0.0;
            }
            // Коэффициенты уравнения производной момента
            a = (q.endValue + r.endValue - q.begValue - r.begValue) * len;
            b = ((q.begValue + r.begValue) * len - m.endValue + m.begValue);
            c = F[0] - m.begValue;
            if ((ABS(a)) < 1e-6) {
                if ((ABS(b)) > 1e-6) {
                    // Может быть еще один экстремум
                    double ks = (-c / b);
                    if (ks > 0.0 && ks < 1.0) {
                        double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
                        double val = -F[1] + F[0] * len * ks -
								len * (m.begValue * ks + (m.endValue - m.begValue) * ks / 2) +
                                ((q.begValue + r.begValue) * ks2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ks3 / 3) * len * len;
						if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
							max = val;
                            ksix = ks;
                        }
                    }
                }
            } else {
                D = b * b - 4 * a * c;
                if (D >= 0.0) {
                    // Могут быть еще 2 экстремума
                    double ks = (-b + sqrt(D)) / (2 * a);
                    if (ks > 0.0 && ks < 1.0) {
                        double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
                        double val = -F[1] + F[0] * len * ks -
                                len * (m.begValue * ks + (m.endValue - m.begValue) * ks / 2) +
                                ((q.begValue + r.begValue) * ks2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ks3 / 3) * len * len;
                        if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
                            max = val;
							ksix = ks;
                        }
                    }
					ks = (-b - sqrt(D)) / (2 * a);
                    if (ks > 0.0 && ks < 1.0) {
                        double ks2 = ks * ks;
                        double ks3 = ks2 * ks;
						double val = -F[1] + F[0] * len * ks -
								len * (m.begValue * ks + (m.endValue - m.begValue) * ks / 2) +
                                ((q.begValue + r.begValue) * ks2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ks3 / 3) * len * len;
                        if (val < min) {
                            min = val;
                            ksin = ks;
                        }
                        if (val > max) {
                            max = val;
                            ksix = ks;
                        }
                    }
                }
            }
            FEMmax.value = max;
            FEMmax.XCoord = EI.X0 + EI.Len * ksix;
            FEMmin.value = min;
			FEMmin.XCoord = EI.X0 + EI.Len * ksin;
}
//---------------------------------------------------------------------------
// Вычисление результатов в заданной точке конечного элемента.
BeamResults BeamFE::resultsAt(double X) {
			double len = EI.Len;
            double ksi = (X - EI.X0) / len;  // Безразмерная координата
            double ksi2 = ksi * ksi;
            double ksi3 = ksi2 * ksi;
			BeamResults br;
            //BeamResults br = new BeamResults();
            // 1. Горизонтальное перемещение
			br.u = dsp[2] + (dsp[5] - dsp[2]) * ksi;
            // 2. Вертикальное перемещение
			br.v = dsp[0] * (1 - 3 * ksi2 + 2 * ksi3) + dsp[1] * len * (ksi - 2 * ksi2 + ksi3) +
                    dsp[3] * (3 * ksi2 - 2 * ksi3) + dsp[4] * len * (ksi3 - ksi2);
            // 3. Угол поворота
            br.c = dsp[0] * (6 * ksi2 - 6 * ksi) / len + dsp[1] * (1 - 4 * ksi + 3 * ksi2) +
                    dsp[3] * (6 * ksi - 6 * ksi2) / len + dsp[4] * (3 * ksi2 - 2 * ksi);
            // 4. Нормальная сила
			br.N = -F[2] - (N.begValue * ksi + (N.endValue - N.begValue) * ksi2 / 2) * len;
            // 5. Поперечная сила
            br.Q = F[0] + ((q.begValue + r.begValue) * ksi + (q.endValue + r.endValue - q.begValue - r.endValue) * ksi2 / 2) * len;
			// 6. Изгибающий момент
			/*
			br.Mz = -F[1] + F[0] * len * ksi + (N.endValue + N.begValue)*(dsp[3]-dsp[0])*0.5*ksi -
					len * (m.begValue * ksi + (m.endValue - m.begValue) * ksi / 2) +
					((q.begValue + r.begValue) * ksi2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ksi3 / 3) * len * len;
			 */
			 /*
			 br.Mz = -F[1] + F[0] * len * ksi -
					len * (m.begValue * ksi + (m.endValue - m.begValue) * ksi / 2) +
					((q.begValue + r.begValue) * ksi2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ksi3 / 3) * len * len;
			 */
			 br.Mz = -F[1] + F[0] * len * ksi  +
					((q.begValue + r.begValue) * ksi2 / 2 + (q.endValue + r.endValue - q.begValue - r.begValue) * ksi3 / 3) * len * len;
			// 7. Реакция упругого основания.
			//br.R = -br.v * kpe.begValue;
            return br;
}
//-------------------------------------------------------------------------------------
// Класс описания данных отдельного региона балки. В пределах региона
// физические параметры меняются линейно.
BeamRegion::BeamRegion() {
          X0 = 0;
          Len = 0;
          begValue = 0;
          endValue = 0;
}

        /**
         * Конструктор региона заданной длины.
		 * @param len
         */
        //void BeamRegion::BeamRegion_a(double len) {
        //    Len = len;
        //}

        /**
         * Конструктор региона, являющегося частью другого региона
         * @param base Базовый регион.
         * @param x0 Координата начала созаваемого региона в ОСК.
         * @param len Длина создаваемого региона.
         */
        //void BeamRegion::BeamRegion_b(BeamRegion *new_region, BeamRegion base, double x0, double len) {
		void BeamFE::BeamRegion_b(BeamRegion *new_region, BeamRegion base, double x0, double len) {
            new_region->X0 = x0;
            new_region->Len = len;
            new_region->begValue = base.valueAt(x0 - base.X0);
            new_region->endValue = base.valueAt(x0 - base.X0 + len);
        }

        double BeamRegion::valueAt(double x) {
            if (Len == 0.0) {
                return begValue;
            }
            if (x < 0.0) {
                return begValue;
            }
            if (x > Len) {
                return endValue;
            }
			return begValue + ((endValue - begValue) / Len) * x;
        }

        Constraint::Constraint() {
          X = 0.0;
          u = 0.0;
          v = 0.0;
          c = 0.0;
        }

        Load::Load() {
          X = 0.0;
          u = 0.0;
          v = 0.0;
          c = 0.0;
        }

    /**
	 * Класс описания эксремального значения
     */
    ExtremalValue::ExtremalValue () {
         value = 0;
         XCoord = 0;
    }
    //String ExtremalValue::toString() {
    //        return String.format("Значение = %g при X = %g", value, XCoord);
    //    }

    /**
     * Класс результатов расчета для заданного сечения балки.
	 * Результаты включают значения сил, моментов, перемещения и углы
     * поворота сечения.
     */
      BeamResults::BeamResults() {
        N = 0;
        Q = 0;
        Mz = 0;
        R = 0;
        u = 0;
        v = 0;
        c = 0;
      }
    /**
     * Конструктор класса по умолчанию. Объект, созданный этим конструктором,
     * произведет расчет балки при максимальной длине конечного элемента равной
     * одной сотой длины балки.
     */
MultiSupportBeam::MultiSupportBeam() {
		maxElementLength = 0.0;
        ConstConstr = 1e30;     // жесткость бесконечно жесткой опоры
        ConstSingular = 1e-15;  // константа, характеризующая вырожденность
        ErrSingular = 1;
        ErrData = 2;
        beamLength = 0.0;

        ConstrForces = NULL;
        fem = NULL;
}
    /**
     * Конструктор класса. Объект, созданный этим конструктором,
	 * произведет расчет балки при максимальной длине конечного элемента заданной
     * в параметре конструктора. Используется, для принудительного управления
     * точностью решения.
     * @param maxElementLength Максимальная длина конечного элемента.
     */

void MultiSupportBeam::MultiSupportBeam_(int N_break0, double maxElemLen, double H_beam0, double Zct_beam0, double I_beam0, double F_beam0) {
		maxElementLength = maxElemLen;
		H_beam = H_beam0;   // Высота сечения
		Zct_beam = Zct_beam0;  // координта ц.т. сечения относительно внешней грани
		I_beam = I_beam0;  // момент инерции сечения
		F_beam = F_beam0;  // площадь сечения
		N_break = N_break0;
}
//------------------------------------------------------------------------------
// Нормализация регионов
int MultiSupportBeam::normaliseRegions(BeamRegion* EI, BeamRegion* N_, BeamRegion* qy,
									   BeamRegion* mz,
									   //BeamRegion* ky,
									   Load* loads, Constraint* constraints, double FE_beam) {
        // Вычисляем длину балки
        beamLength = 0;
        for (int i = 0; i < EI_length; i++) {
			beamLength += EI[i].Len;
        }
        BeamRegion *EF = new BeamRegion[1];
        EF_length = 1;
        EF[0].Len = beamLength;
		EF[0].begValue = FE_beam;
		EF[0].endValue = FE_beam;
		int MaxNormRegion = EI_length + EF_length + q_length + m_length + N_length + kp_length + loads_length + constr_length;
		MaxNormRegion *=10;  // на всякий случай
		// Формируем регионы
		BeamRegion *EIv = new BeamRegion[MaxNormRegion];
        BeamRegion *EFv = new BeamRegion[MaxNormRegion];
        BeamRegion *qxv = new BeamRegion[MaxNormRegion];
        BeamRegion *qyv = new BeamRegion[MaxNormRegion];
		BeamRegion *mzv = new BeamRegion[MaxNormRegion];
		//BeamRegion *kyv = new BeamRegion[MaxNormRegion];
        Load *loadv = new Load[MaxNormRegion];
        Constraint *cnv = new Constraint[MaxNormRegion];
        double x = 0.0;
        int iEI = 0;
        int iEF = 0;
        int iqx = 0;
        int iqy = 0;
		int imz = 0;
		//int iky = 0;
        int ild = 0;
        int icn = 0;

        int kEI = 0;
		int kEF = 0;
        int kqx = 0;
        int kqy = 0;
		int kmz = 0;
		//int kky = 0;
        int kld = 0;
        int kcn = 0;
        double delta = beamLength * 1e-3;
        // Начальные точки нагрузок и закреплений
        Load ldr;
        if ((ABS(loads[ild].X - x)) < delta) {
            ldr.X = x;
			ldr.c = loads[ild].c;
            ldr.u = loads[ild].u;
            ldr.v = loads[ild].v;
			ild++;
        }
        loadv[kld] = ldr;
        kld++;
        Constraint cnr;
        if ((ABS(constraints[icn].X - x)) < delta) {
			cnr.X = x;
            cnr.c = constraints[icn].c;
            cnr.u = constraints[icn].u;
			cnr.v = constraints[icn].v;
            icn++;
        }
        //cnr.u=-1.0;
        cnv[kcn] = cnr;
		kcn++;
		while (iEI < EI_length) {
            double len = EI[iEI].Len;
			while (iqx < N_length) {
                if (N_[iqx].Len > delta) {
                    len = MIN(len, N_[iqx].Len);
                    break;
                }
                iqx++;
            }
            while (iqy < q_length) {
                if (qy[iqy].Len > delta) {
                    len = MIN(len, qy[iqy].Len);
                    break;
                }
                iqy++;
			}
			//
            while (imz < m_length) {
				if (mz[imz].Len > delta) {
                    len = MIN(len, mz[imz].Len);
                    break;
                }
                imz++;
			}
			/*
			while (iky < kp_length) {
                if (ky[iky].Len > delta) {
                    len = MIN(len, ky[iky].Len);
                    break;
                }
				iky++;
			}
			*/
            if (ild < loads_length) {
				len = MIN(len, loads[ild].X - x);
            }
            if (icn < constr_length) {
                len = MIN(len, constraints[icn].X - x);
            }
            if (len<0){
                return ErrData;
            }
            // EF
            iEF = addRegion(EF_length, len, x, EFv, EF, iEF, &kEF);
            // EI
            iEI = addRegion(EI_length, len, x, EIv, EI, iEI, &kEI);
            // qx
            iqx = addRegion(N_length, len, x, qxv, N_, iqx, &kqx);
			// qy
            iqy = addRegion(q_length, len, x, qyv, qy, iqy, &kqy);
            // mz
			imz = addRegion(m_length, len, x, mzv, mz, imz, &kmz);
            // ky
			//iky = addRegion(kp_length, len, x, kyv, ky, iky, &kky);
            // Нагрузки
            ldr.X = x + len;
            ldr.c = 0;
            ldr.u = 0;
            ldr.v = 0;
            if (ild < loads_length) {
				if ((ABS(loads[ild].X - x - len)) < delta) {
					ldr.c = loads[ild].c;
                    ldr.u = loads[ild].u;
                    ldr.v = loads[ild].v;
                    ild++;
                }
            }
            loadv[kld] = ldr;
            kld++;
            // Закрепления
            //cnr = new Constraint();
            cnr.X = x + len;
            cnr.c = 0;
            cnr.u = 0;
            cnr.v = 0;
            if (icn < constr_length) {
                if ((ABS(constraints[icn].X - x - len)) < delta) {
                    cnr.c = constraints[icn].c;
					cnr.u = constraints[icn].u;
                    cnr.v = constraints[icn].v;
                    icn++;
                }
            }
            cnv[kcn] = cnr;
            kcn++;
            x += len;
        }
        EI_BM = new BeamRegion[kEI];
		EF_BM = new BeamRegion[kEI];
		q_BM = new BeamRegion[kEI];
        p_BM = new BeamRegion[kEI];
		m_BM = new BeamRegion[kEI];
		//kp_BM = new BeamRegion[kEI];
        loads_0 = new Load[kEI+1];
        constr_0 = new Constraint[kEI+1];
        EI_BM_length = kEI;
        int i;
        for (i=0; i<EI_BM_length; i++) {
          EI_BM[i] = EIv[i];
          EF_BM[i] = EFv[i];
          q_BM[i] = qyv[i];
          p_BM[i] = qxv[i];
		  m_BM[i] = mzv[i];
		  //kp_BM[i] = kyv[i];
          loads_0[i] = loadv[i];
          constr_0[i] = cnv[i];
        }
        loads_0[i] = loadv[i];
		constr_0[i] = cnv[i];

        delete(EIv);
        delete(EFv);
        delete(qxv);
        delete(qyv);
		delete(mzv);
        //delete(kyv);
		delete(loadv);
        delete(cnv);
		delete(EF);


        return 0;
}
//
//---------------------------------------------------------------------------------------
int MultiSupportBeam::addRegion(int EI_length, double len, double x, BeamRegion *EIv, BeamRegion *EI, int ind, int *kv) {
        BeamRegion EIr;
        EIr.Len = len;
        EIr.X0 = x;
        if (ind < EI_length) {
            EIr.begValue = EI[ind].begValue;
            EIr.endValue = EI[ind].valueAt(len);
            double delta = beamLength * 1e-3;
            if ((ABS(len - EI[ind].Len)) < delta) {
                ind++;
            } else {
                EI[ind].Len -= len;
                EI[ind].begValue = EIr.endValue;
			}
        }
        EIv[*kv] = EIr;
        (*kv)++;
        return ind;
}
//----------------------------------------------------------------------------------------------------
/**
     * Расчет многоопорной балки на упругом основании. <br>
	 * Балка описывается как набор регионов. Регионы располагаются
     * в массивах последовательно, без пропусков. В процессе расчета значения
     * в массивах исходных данных будут изменены.
     * @param EI Массив регионов изгибной жесткости балки.
     * @param N_ Массив регионов эпюры внутренних осевых сил. Если суммарная
     * длина регионов больше суммы длинн регионов EI, то выходящая за пределы
     * балки эпюра игнорируется. Если суммарная
     * длина регионов меньше суммы длинн регионов EI, то на оставшемся участке балки
     * внутренняя осевая сила принимается равной нулю.
     * @param qy Массив регионов распределенной нагрузки вдоль оси Y. Если суммарная
     * длина регионов больше суммы длинн регионов EI, то выходящая за пределы
     * балки нагрузка игнорируется. Если суммарная
     * длина регионов меньше суммы длинн регионов EI, то остаток балки остается
     * ненагруженным.
     * @param mz Массив регионов распределенного изгибающего момента вокруг оси Z.
     * Если суммарная
     * длина регионов больше суммы длинн регионов EI, то выходящая за пределы
     * балки нагрузка игнорируется. Если суммарная
     * длина регионов меньше суммы длинн регионов EI, то остаток балки остается
     * ненагруженным.
     * @param ky Массив регионов упругости основания вдоль оси Y. Оба значения
	 * жесткости на границах региона должны быть одинаковами.
     * Если суммарная
     * длина регионов больше суммы длинн регионов EI, то выходящая за пределы
     * балки нагрузка игнорируется. Если суммарная
     * длина регионов меньше суммы длинн регионов EI, то остаток балки считается
     * свободным. Отрицательная жесткость означает наличие абсолютно жесткого
     * основания на участке.
     * @param loads Сосредоточенные силовые факторы.
	 * @param constraints Опоры.
     * @return Код завершения расчета. <ul>
     * <li> 0 - Успешное завершение расчета.
     * <li> 1 - Вырожденная матрица жесткости конструкции (нулевые жесткости,
     * недостаточно закреплений или сжимающая осевая сила превысила критическую
     * нагрузку потери устойчивости).
     * <li> 2 - Ошибка в исходных данных (нагрузки или закрепления следуют не
     * в порядке возрастания X координаты).
     * </ul>
     */
int MultiSupportBeam::calculate(BeamRegion* EI,
			BeamRegion* N_, BeamRegion* qy,
			BeamRegion* mz,
			//BeamRegion* ky,
			Load* loads, Constraint* constraints, double FE_beam) {
        int rc;                                              
		// Приводим регионы к единому знаменателю
		rc = normaliseRegions(EI, N_, qy,
							  mz,
							  //ky,
							  loads, constraints, FE_beam);
        if (rc != 0) {
           goto free;
        }
        // Получение решения
		rc = Solver();
free:
        delete(EI_BM);
        delete(EF_BM);
        delete(q_BM);
        delete(p_BM);
		delete(m_BM);
		//delete(kp_BM);
        delete(loads_0);
        delete(constr_0);

        return rc;
}
//---------------------------------------------------------------------------
//     * Возвращает результаты расчета балки в заданном сечении.
//     * @param X Координата сечения от левого конца балки.
//     * @return Результаты расчета.
BeamResults MultiSupportBeam::getResults(double X) {
        BeamResults br;
        for (int i = 0; i < fem_length; i++) {
            if (fem[i].EF.X0 <= X && (fem[i].EF.X0 + fem[i].EF.Len) >= X) {
                return fem[i].resultsAt(X);
            }
        }
        return br;
}
//---------------------------------------------------------------------------
int MultiSupportBeam::Solver() {
        int rc;
        double length_break;
		if (maxElementLength == 0.0) {
			length_break = beamLength /100; // размер элемента при разбиении
		} else {
			length_break = maxElementLength; // размер элемента при разбиении
        }

		rc = Solver_Mesh(N_break, length_break);
        return rc;
}
//---------------------------------------------------------------------------
//
int MultiSupportBeam::Solver_Mesh(int N_break0, double length_break) {

        int i, j, k, m, nm, rc, n_break;
        int N_Reg;
        int N_break, Nb;
        double l_elem;
        int nfr;
        double **load_x;
        double **constr_x;
        double *Xb;

		N_Reg = EI_BM_length;

		//double ke[][] = new double[6][6];      // матрица жесткости
		//double load_ke[] = new double[6];      // приведенная нагрузка в узлах элемента
		double **ke;      // матрица жесткости
		double load_ke[6];   // приведенная нагрузка в узлах элемента
		double **k_front; // матрица фронта
		double **g_factor; // верхний треугольник разложенной матрицы
		double *d_beam;    // вектор перемещений модели
		double *l_front;   // нагрузка фронта
		double *d_front;   // перемещения фронта
		double *l_factor;
		int i_break;
		// Максимальный размер рагиона
		double Len_max = 0;
		for (i = 0; i < N_Reg; i++) {
		  Len_max = MAX(Len_max, EI_BM[i].Len);
		}
		length_break = Len_max/N_break0;

		// Подсчет общего числа разбиений

		N_break = 0;
		for (i = 0; i < N_Reg; i++) {
			i_break = Round(EI_BM[i].Len / length_break, 0);
			N_break += MAX(1.0, i_break);
		}

		fem_length = N_break;
		//fem = new BeamFE[N_break];
		if (fem!=NULL) {
		   free(fem);
		}
		fem = (BeamFE*)calloc(N_break, sizeof(BeamFE));

		// Число неизвестных
		nfr = (N_break + 1) * 3;
		// Порядок МЖ - ширина фронта
		nm = 6;

		ke = allocmat(nm, nm);
		load_x = allocmat(N_break + 1, 3);
		constr_x = allocmat(N_break + 1,3);
		k_front = allocmat(nm, nm);
		g_factor = allocmat(nfr, nm);
        //l_beam = new double[nfr];
		Xb = new double[N_break + 1];
        d_beam = new double[nfr];
        l_front = new double[nm];
        d_front = new double[nm / 2];
		l_factor = new double[nfr];


		Xb[0] = 0.0;
		load_x[0][0] = -loads_0[0].v * F1000;
		load_x[0][1] = -loads_0[0].c * F1000*L1000;
		load_x[0][2] = loads_0[0].u * F1000;
		//load_x[0][2] = 0;
		//constr_x[0][0] = constr_0[0].v;
		//constr_x[0][1] = constr_0[0].c;
		//constr_x[0][2] = constr_0[0].u;
		constr_x[0][0] = constr_0[0].v;
		constr_x[0][1] = constr_0[0].c;
		constr_x[0][2] = constr_0[0].u;
		// Разбиение
		Nb = 1;
		for (i = 0; i < N_Reg; i++) {
			i_break = Round(EI_BM[i].Len / length_break, 0);
			n_break =  MAX(1.0, i_break);

			l_elem = EI_BM[i].Len / n_break;
			l_elem *= L1000;  // Преобразование размера элемента в мм
			for (k = 0; k < n_break; k++) {
				Xb[Nb + k] = Xb[Nb + k - 1] + l_elem;
				fem[Nb + k - 1] = BeamFE(Xb[Nb + k - 1], l_elem, EI_BM[i], EF_BM[i], q_BM[i], p_BM[i], m_BM[i]); //, kp_BM[i]);
				load_x[Nb + k][0] = 0.0;
				load_x[Nb + k][1] = 0.0;
				load_x[Nb + k][2] = 0.0;
				constr_x[Nb + k][0] = 0;
				constr_x[Nb + k][1] = 0;
				constr_x[Nb + k][2] = 0;
			}
			// Уточнение координаты конца участка
			Xb[Nb + n_break - 1] = Xb[Nb - 1] + EI_BM[i].Len * L1000;    // длина - в мм
			load_x[Nb + k - 1][0] = -loads_0[i + 1].v * F1000;
			load_x[Nb + k - 1][1] = -loads_0[i + 1].c * F1000*L1000;
			load_x[Nb + k - 1][2] = loads_0[i + 1].u * F1000;
			//constr_x[Nb + k - 1][0] = constr_0[i + 1].v;
			//constr_x[Nb + k - 1][1] = constr_0[i + 1].c;
			//constr_x[Nb + k - 1][2] = constr_0[i + 1].u;
			constr_x[Nb + k - 1][0] = constr_0[i + 1].v;
			constr_x[Nb + k - 1][1] = constr_0[i + 1].c;
			constr_x[Nb + k - 1][2] = constr_0[i + 1].u;
			Nb += n_break;
		}
		// Решение
		// Обнулить матрицу и нагрузки фронта
		for (i = 0; i < nm; i++) {
			l_front[i] = 0.0;
			for (j = 0; j<nm; j++)
			  k_front[i][j] = 0.0;
        }
        // Цикл по числу разбиений
        for (k = 0; k < N_break; k++) {
            //load_ke = fem[k].load;
            // Поместить матрицу во фронт
			//prinv(("Номер элемента "+IntToStr(k)).c_str(), 0, NULL);
            //prinv("load_ke", nm, fem[k].load);
			//---- для печати-----
			//for (i=0; i<6; i++)
			//   for (j=0; j<6; j++)
            //      ke[i][j] = fem[k].ke[i][j];
            //prinm0("ke", nm, nm, ke);
            //------------------------
            m = k * nm / 2;
            for (i = 0; i < nm; i++) {
				for (j = 0; j < nm; j++) {
					k_front[i][j] += fem[k].ke[i][j];
                }
            }
            // подсуммировать к вектору узловых нагрузок узловые нагрузки с элемента
			for (i = 0; i < nm / 2; i++) {
				//load_x[k][i] += fem[k].load[i] * F1000/L1000;
				//load_x[k + 1][i] += fem[k].load[i + nm / 2] * F1000/L1000;
				load_x[k][i] -= fem[k].load[i] * F1000/L1000;
				load_x[k + 1][i] -= fem[k].load[i + nm / 2] * F1000/L1000;
			}
            // подсуммировать к вектору нагрузок фронта узловые нагрузки и наложить связи
			for (i = 0; i < nm / 2; i++) {
                l_front[i] += load_x[k][i];
                if (constr_x[k][i] < 0) {
                    k_front[i][i] += ConstConstr;
                } else {
                    k_front[i][i] += constr_x[k][i];
                }
            }
			//prinm0("k_front перед конденсацией", nm, nm, k_front);
			//prinv("l_front перед конденсацией", nm, l_front);

			rc = mcond_front(nm, nm / 2, k_front, l_front);
            if (rc > 0) {
                //System.out.printf(String.format("\nШаг конденсации %d, матрица вырождена", k+1));
				return rc;
            }
            //prinm0("k_front после конденсациии", nm, nm, k_front);

            // Записать разложенную матрицу и результат прямого хода по нагрузкам
			m = k * nm / 2;
			for (i = 0; i < nm / 2; i++) {
                l_factor[i + m] = l_front[i];
                for (j = i; j < nm; j++) {
					g_factor[i + m][j] = k_front[i][j];
                }
			}
			// Переписать конденсированную матрицу
			m = nm / 2;
			for (i = 0; i < m; i++) {
				l_front[i] = l_front[i + m];
				l_front[i + m] = 0;
				for (j = 0; j < m; j++) {
					k_front[i][j] = k_front[i + m][j + m];
					k_front[i + m][j + m] = 0.0;
					k_front[i][j + m] = 0.0;
					k_front[i + m][j] = 0.0;
				}
			}
			//prinm0("k_front после переписи", nm, nm, k_front);
			//prinv("l_front после переписи", nm, l_front);
		}
		// подсуммировать к вектору нагрузок узловые нагрузки и наложить связи в последнем узле
		for (i = 0; i < 3; i++) {
			l_front[i] += load_x[k][i];
			if (constr_x[k][i] < 0) {
				k_front[i][i] += ConstConstr;
			} else {
				k_front[i][i] += constr_x[k][i];
			}
		}
		// Решение - перемещения в последнем узле
		//prinv("l_front перед решением", nm/2, l_front);
		rc = msolv(nm / 2, k_front, l_front, d_front);
		if (rc > 0) {
			//System.out.printf(String.format("\nРешение, матрица вырождена"));
			return rc;
		}
		//System.out.printf(String.format("\nПеремещения конца последнего участка: vn=%g, fin=%g, un=%g", d_front[0], d_front[1], d_front[2]));
		// Обратный ход - вычисление перемещений балки
		for (i = 0; i < nm / 2; i++) {
			d_beam[i + nfr - nm / 2] = d_front[i];
        }
		mback_front(N_break, nm / 2, g_factor, l_factor, d_beam);
        //prinv("l_beam", nfr, d_beam);

		n_break = (int) MAX(1.0, EI_BM[0].Len / length_break);
        //System.out.printf(String.format("\nПеремещения конца первого участка:    v1=%g, fi1=%g, u1=%g", d_beam[n_break*3],
        //                   d_beam[n_break*3+1], d_beam[n_break*3+2]));
		//-----------------------------------------------------------
		// Формирование реакций в связях
        int n_constr_force = 0;
        // Подсчет количества реакций
        for (i = 0; i < N_break + 1; i++) {
            if (constr_x[i][0] != 0 || constr_x[i][1] != 0 ) {
                n_constr_force++;
            }
        }
		constr_force_length = n_constr_force;

        if (ConstrForces!=NULL) {
           delete(ConstrForces);
        }
        ConstrForces = new Load[n_constr_force];

        for (i = 0, k = 0; i < N_break + 1; i++) {
            if (constr_x[i][0] != 0 || constr_x[i][1] != 0 ) {
                ConstrForces[k].X = Xb[i];
                if (constr_x[i][0] < 0) {
                    ConstrForces[k].v = -d_beam[i * 3] * ConstConstr;
                } else {
                    ConstrForces[k].v = -d_beam[i * 3] * constr_x[i][0];
                }
                if (constr_x[i][1] < 0) {
                    ConstrForces[k].c = -d_beam[i * 3 + 1] * ConstConstr;
                } else {
                    ConstrForces[k].c = -d_beam[i * 3 + 1] * constr_x[i][1];
                }
                //if (constr_x[i][2] < 0) {
				//    ConstrForces[k].u = -d_beam[i * 3 + 2] * ConstConstr;
                //} else {
                //    ConstrForces[k].u = -d_beam[i * 3 + 2] * constr_x[i][2];
                //}
                k++;
            }
        }

		for (i = 0; i < N_break; i++) {
			for (j = 0; j < 6; j++) {
                fem[i].dsp[j] = d_beam[i * 3 + j];
            }
            fem[i].beam_results();

            MaxExtremValue(i, &umax, &(fem[i].FEumax));
            MinExtremValue(i, &umin, &(fem[i].FEumin));
            MaxExtremValue(i, &vmax, &(fem[i].FEvmax));
            MinExtremValue(i, &vmin, &(fem[i].FEvmin));
            MaxExtremValue(i, &cmax, &(fem[i].FEcmax));
            MinExtremValue(i, &cmin, &(fem[i].FEcmin));
            MaxExtremValue(i, &Nmax, &(fem[i].FENmax));
            MinExtremValue(i, &Nmin, &(fem[i].FENmin));
            MaxExtremValue(i, &Qmax, &(fem[i].FEQmax));
            MinExtremValue(i, &Qmin, &(fem[i].FEQmin));
            MaxExtremValue(i, &Mmax, &(fem[i].FEMmax));
            MinExtremValue(i, &Mmin, &(fem[i].FEMmin));
        }
        delete(Xb);
        delete(d_beam);
		delete(l_front);
		delete(d_front);
        delete(l_factor);
        //System.out.printf(String.format("\n"));
        return 0;
}
//-----------------------------------------------------------------------------------------
    /**
	 * Конденсация фронта
	 * @param n0 - порядок исходной матрицы
     * @param m0 - порядок результирующей матрицы
     * @param mat0 - матрица фронта
     * @param load0
     * @return
     */
int MultiSupportBeam::mcond_front(int n0, int m0, double **mat0, double *load0) {
        int i, j, k;
        int l;
        double ss, dd;

        l = n0 - m0;

		/* разложение  */

        for (i = 0; i < l; i++) {
            for (j = 0, ss = 0; j < i; j++) {
                ss += mat0[j][i] * mat0[j][i];
            }
            if (mat0[i][i] - ss < ConstSingular) {
                return ErrSingular;
            }
			dd = mat0[i][i] = sqrt(mat0[i][i] - ss);
            for (k = i + 1; k < l; k++) {
                for (j = 0, ss = 0; j < i; j++) {
                    ss += mat0[j][i] * mat0[j][k];
                }
                mat0[i][k] = (mat0[i][k] - ss) / dd;
            }
		}
        /* прямой ход по матрице k12 */
        for (i = l; i < n0; i++) {
            for (j = 0; j < l; j++) {
                for (k = 0, ss = 0; k < j; k++) {
                    ss += mat0[k][j] * mat0[k][i];
                }
                mat0[j][i] = (mat0[j][i] - ss) / mat0[j][j];
            }
        }
        // прямой ход по нагрузкам
        for (j = 0; j < l; j++) {
            for (k = 0, ss = 0; k < j; k++) {
				ss += mat0[k][j] * load0[k];
            }
            load0[j] = (load0[j] - ss) / mat0[j][j];
        }
        // конденсация матрицы
        for (i = l; i < n0; i++) {
            for (j = l; j < n0; j++) {
                ss = 0;
                for (k = 0; k < l; k++) {
                    ss += mat0[k][i] * mat0[k][j];
				}
                mat0[i][j] -= ss;
            }
        }
        // конденсация нагрузок
        for (j = l; j < n0; j++) {
			ss = 0;
            for (k = 0; k < l; k++) {
                ss += load0[k] * mat0[k][j];
            }
            load0[j] -= ss;
        }
        return 0;
    }  // mcond_front

    /**
     * Обратный ход по нагрузкам
     * @param n_break
     * @param m
	 * @param g_factor
     * @param l_factor
     * @param d_beam
     */
    void MultiSupportBeam::mback_front(int n_break, int m, double **g_factor, double *l_factor, double *d_beam) {
        int i, j, k;
        double ss;

        for (k = m * (n_break - 1); k >= 0; k -= m) {
            for (i = m - 1; i >= 0; i--) {
                ss = 0;
				for (j = i + 1; j < 2 * m; j++) {
                    ss += g_factor[k + i][j] * d_beam[k + j];
				}
                d_beam[k + i] = (l_factor[k + i] - ss) / g_factor[k + i][i];
            }
		}
        return;
    }  // mback_front

    /**
     * треугольное разложение + решение
     * @param nfr -  порядок матрицы
     * @param k1  -  матрица
     * @param rst -  правая часть
     * @param bu  -  вектор неизвестных
     * @return
     */
    int MultiSupportBeam::msolv(int nfr, double **k1, double *rst, double *bu) {
		int i, j, k;
        double ss, dd;

        for (i = 0; i < nfr; i++) {
            ss = 0;
            for (j = 0; j < i; j++) {
                ss += k1[j][i] * k1[j][i];
            }
            ss = k1[i][i] - ss;
            if (ss < ConstSingular) {
                return ErrSingular;
            }
			k1[i][i] = dd = sqrt(ss);
            for (k = i + 1; k < nfr; k++) {
                ss = 0;
                for (j = 0; j < i; j++) {
					ss += k1[j][i] * k1[j][k];
                }
                k1[i][k] = (k1[i][k] - ss) / dd;
            }
        }
        /* прямой ход */
        for (i = 0; i < nfr; i++) {
            ss = 0;
            for (j = 0; j < i; j++) {
                ss += k1[j][i] * bu[j];
            }
            bu[i] = (rst[i] - ss) / k1[i][i];
        }
		for (i = nfr - 1; i >= 0; i--) {
            ss = 0;
            for (j = i + 1; j < nfr; j++) {
                ss += k1[i][j] * bu[j];
            }
            bu[i] = (bu[i] - ss) / k1[i][i];
        }
        return 0;
    }  /* msolv */
//----------------------------------------------------------------------
double **allocmat(int n, int m) {
   double **Xb;
   Xb = new double*[n];
   for (int i=0; i<n; i++) {
     Xb[i] = new double[m];
   }
   return Xb;
}
//---------------------------------------------------------------------------
void MultiSupportBeam::MaxExtremValue(int i, ExtremalValue* max, ExtremalValue *FEmax) {
   if (i==0) {
       max->value = FEmax->value;
       max->XCoord = FEmax->XCoord;
   }
   else {
     if (max->value < FEmax->value) {
         max->value = FEmax->value;
         max->XCoord = FEmax->XCoord;
     }
   }
}
//---------------------------------------------------------------------------
void MultiSupportBeam::MinExtremValue(int i, ExtremalValue* min, ExtremalValue *FEmin) {
   if (i==0) {
      min->value = FEmin->value;
      min->XCoord = FEmin->XCoord;
   }
   else {
     if (min->value > FEmin->value) {
         min->value = FEmin->value;
         min->XCoord = FEmin->XCoord;
     }
   }
}

#pragma package(smart_init)
