Дифференцирование полученного МКЭ решения

[Geckelberryfinn]

Прошу дать совет при решении следующей проблемы.
Есть решение системы дифференциальных уравнение, полученное методом КЭ. Однако, для того, чтобы получить решение исходной задачи, необходимо его дифференцировать. Если пробовать в лоб - т.е. дифференцировать аналитически функции формы и попытаться восстановить решение, получается "не очень" хороший результат (см. рисунки ниже).
Вопрос в слудющем: может быть как-то по-другому осуществлять дифференцирование?

Это решение получено МКЭ


А это результат дифференцирования

[Zai]

У Вас вид решения на первом рисунке достаточно хороший. Нужно проверить Вашу процедуру дифференциирования на одном элементе. Возьмите прямоугольник или прямоугольный треугольник. В нем вручную определите производные а затем сравните их результатами Вашей программы.
Производные в элементах не непрерывные при первом порядке аппроксимации и их можно осреднять( например создать узловую переменную и провести осреднение с весом площади элемента по всем примыкающим к узлу элементам).
В механике деформирования производные от перемещений испльзуются для вычисления напряжений, если они недостаточно непрерывны то нужно дробить сетку.

[VLarin]

Неплохой обзор методов вычисления градиентов (производных) решения в МКЭ есть в этих статьях:
1. Omeragic D., Silvester P.P. Progress in differentiation of approximate data // IEEE Antennas and Propagation Magazine. – 1996. – Vol. 38, N 1, February.
2. Omeragic D., Silvester P.P. Numerical differentiation in magnetic field postprocessing // Int. J. Numer. Model. – 1996, N 9, P. 99 – 113.

Очень часто для вычисления градиентов используют усреднение значений прямого дифференцирования функций формы с весовыми функциями, но ето на мой взгляд не лучший подход:))

Можно воспользоваться Super сonvergence Patch Recovery (Zienkiewicz-Zhu):
Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 1: The recovery technique // Int. j. numer. methods eng. – 1992. – N 33, P. 1331–1364.
Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity // Int. j. numer. methods eng. – 1992. – N 33, P. 1365–1382.

Методов нахождения градиентов куча, а по сему встречный вопрос - как Вы рассчитываете производные решения?

[Geckelberryfinn]

Да, за ссылки спасибо. Посмотрю.
Вычисляю производные решения очень просто. Вот кусок кода

function restoreSolutionDiff(x,y:RealType;XorY: char;RestoreType:byte=0):RealType;
var i,j:integer;
    ra,rb,rc,rd,re,rf:RealType;
begin
Result:=0;
for i:=0 to BlockUsed do
   for j:=0 to High(Blocks[i].Dots) do
      begin
         GetBasisCoeffs(Blocks[i].dots[j],i,ra,rb,rc,rd,re,rf);
         Result:=Result+fmatrix[2*Blocks[i].dots[j]+RestoreType][0]*
            Iif(xory='x', PartXDiff(x,y,ra,rb,rc,rd,re,rf),PartYDiff(x,y,ra,rb,rc,rd,re,rf));
     end;
end;

RestoreType может принимать только 0 или 1 - в зависимости от такого, какую переменную дифференцировать (решалась система УвЧ). PartXDiff, PartYDiff - непосредственно функции вычисления частных производных, там предельно все просто:

f
unction PartXDiff(x,y,a,b,c,d,e,f:realtype):RealType;
begin
Result:=b+d*y+2*e*x;
end;

function PartYDiff(x,y,a,b,c,d,e,f:realtype):RealType;
begin
Result:=(c+d*x+2*f*y);
end;

т.к. функция формы - $a+b\cdot x+c\cdot y+d\cdot x\cdot y+e\cdot x^2+f\cdot y^2$

[Geckelberryfinn]

Нашел кое-что про SPRecovery в Zienkiewicz O.C., Taylor R.L. Vol. 1. The finite element method. The basis. Но не совсем понял суть этого метода. Там вроде как находятся коэффициенты полинома, аппроксимирующего производную, но мне не совсем понятны данные, на которых эта аппроксимация осуществляется. Никто не сталкивался с таким методом и не мог бы разъяснить его более подробно?

[Zai]

Ничего более продвинутого более чем тероия вычислений среднего значения производной на геометрической фигуре конечного элемента Вы там не сможите найдти. Когда различие производных в решении отличается на несколько десятых процента, самый верный путь погрузиться в чтении статей на несклько месяцев. Если же отличие ничего не имеет с качественным решением, то естественно провести некоторое время в функциональном наполнении программ debugger.
Если решение глобального вектора смещений у Вас физично, то нужно разобраться почему нарушается аппроксимация производных. Это скорее задача поиска ошибок в программе.

[VLarin]

Geckelberryfinn писал(а):

Нашел кое-что про SPRecovery в Zienkiewicz O.C., Taylor R.L. Vol. 1. The finite element method. The basis. Но не совсем понял суть этого метода. Там вроде как находятся коэффициенты полинома, аппроксимирующего производную, но мне не совсем понятны данные, на которых эта аппроксимация осуществляется. Никто не сталкивался с таким методом и не мог бы разъяснить его более подробно?

Идея там следующая. Показывается, что для полиномов Лагранжа значения производных, полученных путем их прямого дифференцирования, совпадают с точными значениями производных в точках на конечных элементе, называемых суперсходимыми. И для полиномов Лагранжа эти точки совпадают с точками Гауссовых квадратур, используемых для интегрирования на КЭ.

Используя значения производных в этих точках можно попробовать восстановить картину производных в окрестности этих точек.
Для этого в SPR используют метод наименьших квадратов - с его помощью определяются значения коэффициентов полиномов, наиболее точно описывающих производные в некоторой рассматриваемой области. После того, как были найдены значения коэффициентов, можно с их помощью определить значения производных в узлах рассматриваемы КЭ.

Процесс восстановления значений производных выглядит следующим образом:
1. Из списка узлов выбирается некоторый узел.
2. Определяются КЭ, окружающие этот узел и принадлежащие одной среде.
3. Для найденных КЭ определяются значения производных путем прямого дифференцирования функций формы в Гауссовых точках.
4. Используя МНК определяются коэффициенты полиномов, описывающих изменение этих производных в окрестности рассматриваемого узла.
5. По найденным полиномам определяются значения востановленные производных в рассматриваемом узле, а также промежуточных узлах рассматриваемы КЭ.
6. Повторяются пп.1-5 для всех остальных не рассмотренных сред (подобластей) и узлов.

После этого получаем значения восстановаленных производных во всех узлах и для всех подобластей. Определение значений производных внутри КЭ производится аналогично
определению значения потенциала (или значения искомой функции, я просто электромагнитными полями занимаюсь:)) - т.е. используя функции формы КЭ.