41double L2Norm(GeometricField<scalar, fvPatchField, volMesh>& field)
44 a = Foam::sqrt(fvc::domainIntegrate(field * field).value());
49double L2Norm(GeometricField<vector, fvPatchField, volMesh>& field)
52 a = Foam::sqrt(fvc::domainIntegrate(field & field).value());
57double LinfNorm(GeometricField<scalar, fvPatchField, volMesh>& field)
60 a = Foam::max(Foam::sqrt(field.internalField() *
61 field.internalField())).value();
66double LinfNorm(GeometricField<vector, fvPatchField, volMesh>& field)
69 Info <<
"LinfNorm(GeometricField<vector, fvPatchField, volMesh>& field) is still to be implemented"
77template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
78double errorFrobRel(GeometricField<Type, PatchField, GeoMesh>& field1,
79 GeometricField<Type, PatchField, GeoMesh>& field2, List<label>* labels)
82 autoPtr<GeometricField<Type, PatchField, GeoMesh>> errField;
83 autoPtr<GeometricField<Type, PatchField, GeoMesh>> field1_S;
84 autoPtr<GeometricField<Type, PatchField, GeoMesh>> field2_S;
85 autoPtr<fvMeshSubset> submesh;
89 submesh = autoPtr<fvMeshSubset>(
new fvMeshSubset(field1.mesh()));
91 submesh->setCellSubset(*labels);
93 submesh->setLargeCellSubset(*labels);
95 GeometricField<Type, PatchField, GeoMesh> field1tmp(submesh->interpolate(
97 GeometricField<Type, PatchField, GeoMesh> field2tmp(submesh->interpolate(
99 field1_S = autoPtr<GeometricField<Type, PatchField, GeoMesh>>
100 (
new GeometricField<Type, PatchField, GeoMesh>(field1tmp.clone()));
101 field2_S = autoPtr<GeometricField<Type, PatchField, GeoMesh>>
102 (
new GeometricField<Type, PatchField, GeoMesh>(field2tmp.clone()));
106 field1_S = autoPtr<GeometricField<Type, PatchField, GeoMesh>>
107 (
new GeometricField<Type, PatchField, GeoMesh>(field1));
108 field2_S = autoPtr<GeometricField<Type, PatchField, GeoMesh>>
109 (
new GeometricField<Type, PatchField, GeoMesh>(field2));
112 errField = autoPtr<GeometricField<Type, PatchField, GeoMesh>>
113 (
new GeometricField<Type, PatchField, GeoMesh>(field1_S() - field2_S()));
128template double errorFrobRel(GeometricField<scalar, fvPatchField, volMesh>&
130 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
131template double errorFrobRel(GeometricField<vector, fvPatchField, volMesh>&
133 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
135template double errorFrobRel(GeometricField<scalar, fvsPatchField, surfaceMesh>&
137 GeometricField<scalar, fvsPatchField, surfaceMesh>& field2,
138 List<label>* labels);
143 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
146 autoPtr<GeometricField<T, fvPatchField, volMesh>> errField;
147 autoPtr<GeometricField<T, fvPatchField, volMesh>> field1_S;
148 autoPtr<GeometricField<T, fvPatchField, volMesh>> field2_S;
149 autoPtr<fvMeshSubset> submesh;
153 submesh = autoPtr<fvMeshSubset>(
new fvMeshSubset(field1.mesh()));
155 submesh->setCellSubset(*labels);
157 submesh->setLargeCellSubset(*labels);
159 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
161 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
163 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
164 (
new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
165 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
166 (
new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
170 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
171 (
new GeometricField<T, fvPatchField, volMesh>(field1));
172 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
173 (
new GeometricField<T, fvPatchField, volMesh>(field2));
176 errField = autoPtr<GeometricField<T, fvPatchField, volMesh>>
177 (
new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
192template double errorLinfRel(GeometricField<scalar, fvPatchField, volMesh>&
194 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
195template double errorLinfRel(GeometricField<vector, fvPatchField, volMesh>&
197 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
200double errorL2Abs(GeometricField<vector, fvPatchField, volMesh>& field1,
201 GeometricField<vector, fvPatchField, volMesh>& field2, volScalarField& Volumes)
204 volScalarField diffFields2 = (((field1 - field2) & (field1 - field2)) *
206 double err = Foam::sqrt(gSum(diffFields2));
211double errorL2Abs(GeometricField<scalar, fvPatchField, volMesh>& field1,
212 GeometricField<scalar, fvPatchField, volMesh>& field2, volScalarField& Volumes)
215 volScalarField diffFields2 = (((field1 - field2) * (field1 - field2)) *
217 double err = Foam::sqrt(gSum(diffFields2));
222double errorL2Abs(GeometricField<T, fvPatchField, volMesh>& field1,
223 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
225 autoPtr<GeometricField<T, fvPatchField, volMesh>> errField;
226 autoPtr<GeometricField<T, fvPatchField, volMesh>> field1_S;
227 autoPtr<GeometricField<T, fvPatchField, volMesh>> field2_S;
228 autoPtr<fvMeshSubset> submesh;
232 submesh = autoPtr<fvMeshSubset>(
new fvMeshSubset(field1.mesh()));
234 submesh->setCellSubset(*labels);
236 submesh->setLargeCellSubset(*labels);
238 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
240 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
242 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
243 (
new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
244 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
245 (
new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
249 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
250 (
new GeometricField<T, fvPatchField, volMesh>(field1));
251 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
252 (
new GeometricField<T, fvPatchField, volMesh>(field2));
255 errField = autoPtr<GeometricField<T, fvPatchField, volMesh>>
256 (
new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
257 double err =
L2Norm(errField());
261template double errorL2Abs(GeometricField<scalar, fvPatchField, volMesh>&
263 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
264template double errorL2Abs(GeometricField<vector, fvPatchField, volMesh>&
266 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
269template<
class T,
template<
class>
class PatchField,
class GeoMesh>
270Eigen::MatrixXd
errorFrobRel(PtrList<GeometricField<T, PatchField, GeoMesh>>&
272 PtrList<GeometricField<T, PatchField, GeoMesh>>& fields2, List<label>* labels)
276 if (fields1.size() != fields2.size())
278 Info <<
"The two fields do not have the same size, code will abort" << endl;
282 err.resize(fields1.size(), 1);
284 for (label k = 0; k < fields1.size(); k++)
286 err(k, 0) =
errorFrobRel(fields1[k], fields2[k], labels);
287 Info <<
" Error is " << err[k] << endl;
294 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields1,
295 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields2,
296 List<label>* labels);
298 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields1,
299 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields2,
300 List<label>* labels);
302 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& fields1,
303 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& fields2,
304 List<label>* labels);
309 PtrList<GeometricField<T, fvPatchField, volMesh>>& fields1,
310 PtrList<GeometricField<T, fvPatchField, volMesh>>& fields2,
311 PtrList<volScalarField>& Volumes)
313 M_Assert(fields1.size() == fields2.size(),
314 "The two fields do not have the same size, code will abort");
315 M_Assert(fields1.size() == Volumes.size(),
316 "The volumes field and the two solution fields do not have the same size, code will abort");
318 err.resize(fields1.size(), 1);
320 for (label k = 0; k < fields1.size(); k++)
322 err(k, 0) =
errorL2Abs(fields1[k], fields2[k], Volumes[k]);
323 Info <<
" Error is " << err[k] << endl;
330 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields1,
331 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields2,
332 PtrList<volScalarField>& Volumes);
334 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields1,
335 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields2,
336 PtrList<volScalarField>& Volumes);
339Eigen::MatrixXd
errorL2Abs(PtrList<GeometricField<T, fvPatchField, volMesh>>&
341 PtrList<GeometricField<T, fvPatchField, volMesh>>& fields2, List<label>* labels)
345 if (fields1.size() != fields2.size())
347 Info <<
"The two fields do not have the same size, code will abort" << endl;
351 err.resize(fields1.size(), 1);
353 for (label k = 0; k < fields1.size(); k++)
355 err(k, 0) =
errorL2Abs(fields1[k], fields2[k], labels);
356 Info <<
" Error is " << err[k] << endl;
363 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields1,
364 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields2,
365 List<label>* labels);
367 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields1,
368 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields2,
369 List<label>* labels);
372double errorL2Rel(GeometricField<T, fvPatchField, volMesh>& field1,
373 GeometricField<T, fvPatchField, volMesh>& field2, List<label>* labels)
376 autoPtr<GeometricField<T, fvPatchField, volMesh>> errField;
377 autoPtr<GeometricField<T, fvPatchField, volMesh>> field1_S;
378 autoPtr<GeometricField<T, fvPatchField, volMesh>> field2_S;
379 autoPtr<fvMeshSubset> submesh;
383 submesh = autoPtr<fvMeshSubset>(
new fvMeshSubset(field1.mesh()));
385 submesh->setCellSubset(*labels);
387 submesh->setLargeCellSubset(*labels);
389 GeometricField<T, fvPatchField, volMesh> field1tmp(submesh->interpolate(
391 GeometricField<T, fvPatchField, volMesh> field2tmp(submesh->interpolate(
393 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
394 (
new GeometricField<T, fvPatchField, volMesh>(field1tmp.clone()));
395 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
396 (
new GeometricField<T, fvPatchField, volMesh>(field2tmp.clone()));
400 field1_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
401 (
new GeometricField<T, fvPatchField, volMesh>(field1));
402 field2_S = autoPtr<GeometricField<T, fvPatchField, volMesh>>
403 (
new GeometricField<T, fvPatchField, volMesh>(field2));
406 errField = autoPtr<GeometricField<T, fvPatchField, volMesh>>
407 (
new GeometricField<T, fvPatchField, volMesh>(field1_S() - field2_S()));
409 if (
L2Norm(field1) <= 1e-6)
422template double errorL2Rel(GeometricField<scalar, fvPatchField, volMesh>&
424 GeometricField<scalar, fvPatchField, volMesh>& field2, List<label>* labels);
425template double errorL2Rel(GeometricField<vector, fvPatchField, volMesh>&
427 GeometricField<vector, fvPatchField, volMesh>& field2, List<label>* labels);
430Eigen::MatrixXd
errorL2Rel(PtrList<GeometricField<T, fvPatchField, volMesh>>&
431 fields1, PtrList<GeometricField<T, fvPatchField, volMesh>>& fields2,
436 if (fields1.size() != fields2.size())
438 Info <<
"The two fields do not have the same size, code will abort" << endl;
442 err.resize(fields1.size(), 1);
444 for (label k = 0; k < fields1.size(); k++)
446 err(k, 0) =
errorL2Rel(fields1[k], fields2[k], labels);
447 Info <<
" Error is " << err[k] << endl;
454 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields1,
455 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& fields2,
456 List<label>* labels);
458 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields1,
459 PtrList<GeometricField<vector, fvPatchField, volMesh>>& fields2,
460 List<label>* labels);
463double H1Seminorm(GeometricField<scalar, fvPatchField, volMesh>& field)
466 a = Foam::sqrt(fvc::domainIntegrate(fvc::grad(field) & fvc::grad(
472double H1Seminorm(GeometricField<vector, fvPatchField, volMesh>& field)
475 a = Foam::sqrt(fvc::domainIntegrate(fvc::grad(field)
476 && fvc::grad(field)).value());
481template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
482double frobNorm(GeometricField<Type, PatchField, GeoMesh>& field)
490template double frobNorm(GeometricField<scalar, fvPatchField, volMesh>& field);
491template double frobNorm(GeometricField<vector, fvPatchField, volMesh>& field);
498 label patchID =
mesh.boundaryMesh().findPatchID(patch);
499 const polyPatch& cPatch =
mesh.boundaryMesh()[patchID];
501 const labelUList& faceCells = cPatch.faceCells();
505 label faceOwner = faceCells[faceI] ;
506 scalar faceArea =
mesh.magSf().boundaryField()[patchID][faceI];
507 L2 += faceArea * field[faceOwner] * field[faceOwner];
509 return Foam::sqrt(L2);
513 List<scalar>& field2, word patch)
515 M_Assert(field1.size() == field2.size(),
516 "The two fields do not have the same size, code will abort");
519 label patchID =
mesh.boundaryMesh().findPatchID(patch);
520 const polyPatch& cPatch =
mesh.boundaryMesh()[patchID];
521 M_Assert(field1.size() == cPatch.size(),
522 "The filds must have the same size of the patch, code will abort");
525 scalar faceArea =
mesh.magSf().boundaryField()[patchID][faceI];
526 L2 += faceArea * field1[faceI] * field2[faceI];
536 label patchID =
mesh.boundaryMesh().findPatchID(patch);
537 const polyPatch& cPatch =
mesh.boundaryMesh()[patchID];
539 const labelUList& faceCells = cPatch.faceCells();
542 label faceOwner = faceCells[faceI] ;
546 Linf = std::abs(field[faceOwner]);
548 else if (std::abs(field[faceOwner]) > Linf)
550 Linf = std::abs(field[faceOwner]);
561 label patchID =
mesh.boundaryMesh().findPatchID(patch);
562 const polyPatch& cPatch =
mesh.boundaryMesh()[patchID];
564 const labelUList& faceCells = cPatch.faceCells();
568 label faceOwner = faceCells[faceI] ;
569 scalar faceArea =
mesh.magSf().boundaryField()[patchID][faceI];
570 integral += faceArea * field[faceOwner];
580 label patchID =
mesh.boundaryMesh().findPatchID(patch);
581 const polyPatch& cPatch =
mesh.boundaryMesh()[patchID];
583 const labelUList& faceCells = cPatch.faceCells();
584 int patchSize =
mesh.magSf().boundaryField()[patchID].size();
585 std::string message =
"The input list (size = " + std::to_string(field.size())
586 +
") must have the same size of the patch mesh (size = "
587 + std::to_string(patchSize) +
")";
588 M_Assert( patchSize == field.size(), message.c_str());
592 label faceOwner = faceCells[faceI] ;
593 scalar faceArea =
mesh.magSf().boundaryField()[patchID][faceI];
594 integral += faceArea * field[faceI];
forAll(example_CG.gList, solutionI)
#define M_Assert(Expr, Msg)
Header file of the ITHACAerror file.
static Eigen::VectorXd field2Eigen(GeometricField< tensor, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
Namespace to implement some useful assign operation of OF fields.
double LinfNorm(GeometricField< scalar, fvPatchField, volMesh > &field)
double errorLinfRel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in Linf norm.
double LinfNormOnPatch(fvMesh &mesh, volScalarField &field, word patch)
Evaluate the Linf norm of a field on a boundary patch.
double H1Seminorm(GeometricField< scalar, fvPatchField, volMesh > &field)
double L2Norm(GeometricField< scalar, fvPatchField, volMesh > &field)
double integralOnPatch(fvMesh &mesh, volScalarField &field, word patch)
Evaluate the integral on a patch.
double L2normOnPatch(fvMesh &mesh, volScalarField &field, word patch)
Evaluate the L2 norm of a field on a boundary patch.
double errorL2Rel(GeometricField< T, fvPatchField, volMesh > &field1, GeometricField< T, fvPatchField, volMesh > &field2, List< label > *labels)
Computes the relative error between two geometric Fields in L2 norm.
double frobNorm(GeometricField< Type, PatchField, GeoMesh > &field)
Evaluate the Frobenius norm of a field.
double errorFrobRel(GeometricField< Type, PatchField, GeoMesh > &field1, GeometricField< Type, PatchField, GeoMesh > &field2, List< label > *labels)
Computes the relative error between two Fields in the Frobenius norm.
double L2productOnPatch(fvMesh &mesh, List< scalar > &field1, List< scalar > &field2, word patch)
Evaluate the L2 inner product between two scalarLists.
double errorL2Abs(GeometricField< vector, fvPatchField, volMesh > &field1, GeometricField< vector, fvPatchField, volMesh > &field2, volScalarField &Volumes)