1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file unit_NedelecEdgeTetrahedron.cpp
19 \author E. Lunéville
20 \since 16 sep 2015
21 \date 16 sep 2015
22
23 Low level tests of Nedelec's edge element on tetrahedron
24 and a global test
25
26 This function may either creates a reference file storing the results (check=false)
27 or compares results to those stored in the reference file (check=true)
28 It returns reporting information in a string
29 */
30
31 #include "xlife++-libs.h"
32 #include "testUtils.hpp"
33
34 #include <iostream>
35 #include <sstream>
36
37 using namespace xlifepp;
38
39 namespace unit_NedelecEdgeTetrahedron
40 {
41
42 typedef GeomDomain& Domain;
43
44 Real omg=1, eps=1, mu=1, a=pi_, ome=omg* omg* mu* eps;
45
f(const Point & P,Parameters & pa=defaultParameters)46 Vector<Real> f(const Point& P, Parameters& pa = defaultParameters)
47 {
48 Real x=P(1), y=P(2), z=P(3);
49 Vector<Real> res(3);
50 Real c=3*a*a-ome;
51 res(1)=-2*cos(a*x)*sin(a*y)*sin(a*z)*c;
52 res(2)= sin(a*x)*cos(a*y)*sin(a*z)*c;
53 res(3)= sin(a*x)*sin(a*y)*cos(a*z)*c;
54 return res;
55 }
56
solex(const Point & P,Parameters & pa=defaultParameters)57 Vector<Real> solex(const Point& P, Parameters& pa = defaultParameters)
58 {
59 Real x=P(1), y=P(2), z=P(3);
60 Vector<Real> res(3);
61 res(1)=-2*cos(a*x)*sin(a*y)*sin(a*z);
62 res(2)= sin(a*x)*cos(a*y)*sin(a*z);
63 res(3)= sin(a*x)*sin(a*y)*cos(a*z);
64 return res;
65 }
66
unit_NedelecEdgeTetrahedron(bool check)67 String unit_NedelecEdgeTetrahedron(bool check)
68 {
69 String rootname = "unit_NedelecEdgeTetrahedron";
70 trace_p->push(rootname);
71 std::stringstream out;
72 out.precision(testPrec);
73 verboseLevel(110);
74
75 Point x(1./3,1./3,1./3);
76 Interpolation* interp;
77 interp=findInterpolation(Nedelec, _firstFamily, 1, Hrot);
78 out<<"\n\n================ Nedelec edge tetrahedron, first family, order 1 (NE1_1) ======================"<<eol;
79 NedelecEdgeFirstTetrahedronPk NE11(interp);
80 NE11.print(out);
81 NE11.computeShapeValues(x.begin(),true);
82 NE11.printShapeValues(out,x.begin());
83 out<<"\n\n================ Nedelec edge tetrahedron, first family, order 2 (NE1_2) ======================"<<eol;
84 interp=findInterpolation(Nedelec, _firstFamily, 2, Hrot);
85 NedelecEdgeFirstTetrahedronPk NE12(interp);
86 NE12.print(out);
87 NE12.computeShapeValues(x.begin(),true);
88 NE12.printShapeValues(out,x.begin());
89 out<<"\n\n================ Nedelec edge tetrahedron, first family, order 3 (NE1_3) ======================"<<eol;
90 interp=findInterpolation(Nedelec, _firstFamily, 3, Hrot);
91 NedelecEdgeFirstTetrahedronPk NE13(interp);
92 NE13.print(out);
93 NE13.computeShapeValues(x.begin(),true);
94 NE13.printShapeValues(out,x.begin());
95 out<<eol;
96 out<<"\n\n================ Nedelec edge tetrahedron, first family, order 4 (NE1_4) ======================"<<eol;
97 interp=findInterpolation(Nedelec, _firstFamily, 4, Hrot);
98 NedelecEdgeFirstTetrahedronPk NE14(interp);
99 NE14.print(out);
100 NE14.computeShapeValues(x.begin(),true);
101 NE14.printShapeValues(out,x.begin());
102 out<<"\n================================================================================================"<<eol<<eol;
103
104 /* test : solve Maxwell harmonic equation on unit cube Omega=]0,1[x]0,1[x]0,1[
105 curl curl(E) -w^2.mu.eps E = -i.w.mu.J=fJ in Omega
106 E x n = 0 on dOmega
107 Solution (a=k.pi) :
108 Ex = -2cos(ax)sin(ay)sin(az)
109 Ey = cos(ay)sin(ax)sin(az)
110 Ez = cos(az)sin(ax)sin(ay)
111 -> fJ= (3a^2-w^2.mu.eps)E
112 */
113
114 //numberOfThreads(1);
115 Number ordmin=1, ordmax=1;
116 Number na=10;
117
118 for(Number ord=ordmin; ord<=ordmax; ord++)
119 {
120 Number nh=1+(na/ord);
121 Mesh mesh3d(Cube(_origin=Point(0.,0.,0.),_length=1, _nnodes=nh, _domain_name="Omega",_side_names=Strings(6,"Gamma")),_tetrahedron,1,_gmsh);
122 Domain omega=mesh3d.domain("Omega");
123 Domain gamma=mesh3d.domain("Gamma");
124
125 out<<"=============== test Nedelec edge first family of order "<<ord<<" ==============="<<eol;
126 //create spaces
127 Space V(omega,interpolation(_Nedelec,_firstFamily,ord,Hrot),"V",false);
128 Unknown e(V,"E"); TestFunction q(e,"q");
129 Space W(omega,interpolation(_Lagrange,_standard,ord,H1),"W",true);
130 Unknown u(W,"u",3); TestFunction v(u,"v"); //for H1 projection
131 out<<"number of elements on edge : "<<nh-1;
132 out<<" number of tetrahedra : "<<mesh3d.nbOfElements();
133 out<<" number of dofs : "<<V.nbDofs()<<eol;
134 elapsedTime("mesh and space", std::cout);
135
136 //create FE terms
137 BilinearForm aev=intg(omega,curl(e)|curl(q))-ome*intg(omega,e|q);
138 LinearForm l=intg(omega,f|q);
139 EssentialConditions ecs = (ncross(e) | gamma)=0;
140 TermMatrix A(aev, ecs,"A");
141 TermVector b(l,"B");
142 elapsedTime("FE computation", std::cout);
143
144 //solve system
145 TermVector E=directSolve(A,b);
146 elapsedTime("solving", std::cout);
147
148 //interpolate "edge" solution
149 TermVector Ei=interpolate(u,omega,E,"Ei");
150 saveToFile("Ei_"+tostring(ord),Ei,_vtu);
151
152 // Pk interpolation, L2 projection
153 TermMatrix N(intg(omega,e|v),"N");
154 TermMatrix M(intg(omega,u|v),"M");
155 TermVector NE=N*E;
156 TermVector EP1=directSolve(M,NE,true); //L2 projection
157
158 //compute error
159 TermVector E_ex(u,omega,solex);
160 TermVector E_er=E_ex-EP1;
161 out<<"L2 error on Lagrange projection= "<<sqrt(real((M*E_er|E_er)))<<eol<<eol;
162 elapsedTime("error computation", std::cout);
163 //saveToFile("Ep_"+tostring(ord),EP1,_vtu);
164 //saveToFile("E_ex",E_ex,_vtu);
165 // TermVector Eg=E|gamma;
166 // theCout<<Eg<<eol;
167 }
168
169 //------------------------------------------------------------------------------------
170 // save results in a file or compare results with some references value in a file
171 //------------------------------------------------------------------------------------
172 trace_p->pop();
173 if(check) { return diffResFile(out, rootname); }
174 else { return saveResToFile(out, rootname); }
175 }
176
177 }
178